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I.  Program  Overview 

A  memory  is  any  physical  system  with  a  large  number  of  distinguishable  states.  Through  the  work 
undertaken  in  this  program  we  have  shown  that  an  optical  field  can  be  used  to  distinguish  among  N 
different  physical  configurations  of  a  nanostructure  whose  volume  is  comparable  to  a  cubic  wavelength, 
thus  achieving  log2(N)  bits  of  capacity  within  a  single  resolution  element.  We  have  studied  the  use  of  near- 
field  optical  measurement  together  with  spatially-  and  spectrally-selective  defects  to  define  and  readout 
sub-resolution  nanostructure  storage  configurations.  We  have  successfully  identified,  studied,  and 
characterized  nanostructure  configurations  that  provide  optically  distinguishable  states  with  large  interstate 
distances.  During  year  1  we  focused  on  sub-resolution  surface-relief  structures  combined  with  near-field 
detector  arrays  to  demonstrate  approximate  .storage  densities  of  25  bits/pm2.  During  year  2  we  extended 
this  work  to  so-called  nano-structured  voxels  in  which  we  employ  sub-resolution  volumetric  degrees  of 
freedom  and  near-field  detector  arrays.  We  design  the  data  carrying  volumetric  defects  so  as  to  exploit  a 
plasmon  resonance  in  their  metallic  nanoshells  and  we  predict  significantly  improved  storage  densities  of 
nearly  300  bits/pm2.  This  work  has  resulted  in  an  improved  understanding  of  the  interaction  between  space- 
time  electromagnetic  fields  and  various  (imperfect)  nanostructure  volumes,  new  near-field  optical 
characterization  methods  to  determine  the  sub-resolution  3D  configuration  of  artificial  nanostructures,  and 
new  bounds  on  the  abilities  of  optical  fields  to  probe  physical  processes  on  length  scales  below  the  optical 
wavelength. 

One  unanticipated  outcome  of  the  work  outlined  above  was  the  creation  of  a  novel  method  of  designing 
nonlinear  codes.  These  codes  were  necessary  in  order  to  tolerate  the  crosstalk  that  characterizes  nano- 
structured  optical  memories.  We  invented  these  codes  for  use  in  the  memory  application  and  later  found 
them  to  be  valuable  for  use  in  fiber  communications. 

Our  research  activities  within  this  program  have  progressed  along  two  distinct  lines  of  inquiry.  The  original 
program  vision  was  exclusively  concerned  with  static  optical  memory  structures.  Within  this  component  of 
the  program  we  have  analyzed  and  demonstrated  the  utilization  of  sub-resolution  degrees  of  freedom  to 
achieve  room-temperature  high-density  optical  storage  as  outlined  above.  A  complementary  component  of 
our  activity  grew  out  of  some  early  success  and  gave  rise  to  a  parallel  research  path.  This  parallel  path  has 
been  concerned  with  the  use  of  sub-resolution  volumetric  degrees  of  freedom  to  control  light  in  the  lime 
domain.  This  form  of  control  facilitates  simple  delay-based  storage  and  can  be  important  in  pulse  shaping 
and  optical  buffering  applications.  During  year  2  we  pursued  a  theoretical  understanding  of  information 
propagation  in  such  systems.  We  also  experimentally  demonstrated  that  the  information  velocity  in  so- 
called  fast-light  media  remains  below  the  limit  defined  by  relativistic  causality.  Our  year  3  activities 
continued  along  these  lines.  We  extended  our  theoretical  understanding  of  information  velocity  into  the 
domain  of  slow-light  media  and  we  have  experimentally  demonstrated  the  surprising  result  that  even  for 
systems  in  which  the  group  velocity  is  much  less  than  the  speed  of  light  in  vacuum  the  information  velocity 
remains  close  to  c.  The  no-cost-extension  of  this  work  has  concluded  with  a  novel  method  of  compensating 
distortion  in  slow-light  optical  delay  lines. 

This  report  is  organized  as  follows.  Our  research  is  conveniently  described  in  terms  of  8  different  projects. 

A  brief  summary  of  each  project  and  its  outcomes  is  provided  below.  For  work  that  has  already  resulted  in 
one  or  more  peer-reviewed  journal  articles  these  are  included  as  a  detailed  description  of  the  relevant 
technical  details. 


II.  Surface-Relief  Nano-Structures 

Our  main  project  component  has  dealt  with  the  study  of  surface  relief  structures.  Within  this  component 
we  have  made  contributions  in  three  areas:  analysis  and  simulation,  planar  memory  concepts,  and  detector 
designs  and  realizations.  We  continued  and  extended  the  development  of  our  FDTD  simulators  to  study  the 
memory  density  of  nanoscale  surface-relief  structures.  We  have  advanced  several  approaches  to  calculate 
the  clique-based  capacity  of  these  structures.  We  have  considered  and  bounded  statistical  cross-talk 
between  different  relief  configurations.  We  have  also  considered  a  number  of  new  planar  memory 
concepts.  These  have  included  the  study  of  relief  layers  consisting  of  polymer,  a  low  (n=l.8)  contrast 
medium,  versus  semiconductor,  a  high  (nK3.5)  contrast  medium;  structures  that  represent  binary  versus 
non-binary  states;  and  relief  structures  that  emphasize  hole  versus  backfill  configurations.  Parameterization 
of  a  variety  of  the  relief  parameters  including  the  gap  size  between  the  cells  of  interest,  the  base  height,  and 
the  relief  heights  has  been  accomplished  to  determine  their  influence  on  the  number  of  distinguishable 
states  and,  hence,  the  capacity.  More  details  about  this  work  can  be  found  in  Appendix  1 . 

III.  Nano-Sixuctured  Voxels 

Our  second  main  project  component  has  dealt  with  the  study  of  volumetric  memory  structures.  Within  this 
component  we  have  made  contributions  in  two  main  areas:  analysis  and  simulation  of  nano-structured 
voxels  (NSVs)  defined  with  dielectric,  metallic,  and/or  plasmonic  elements  and  the  use  of  optical  tweezers 
for  novel  fabrication  methods  to  realize  these  volumetric  structures.  We  have  used  our  FDTD  simulator  to 
find  which  nano-structured  configurations  can  be  used  to  realize  high-density  memory  storage  devices. 
Detector  arrays  were  placed  appropriately  under  the  scattering  elements.  Very  distinctive  states  have  been 
realized  with  the  metallic  and  plasmonic  scattering  elements.  We  have  found  that  the  dielectric,  metallic, 
and/or  plasmonic  elements  could  be  positioned  in  the  requisite  matrix  using  optical  tweezer  technologies. 
This  fabrication  approach  is  being  analyzed  further  to  understand  its  limitations  and  practicality. 

IV.  Nonlinear  Coding  for  Fiber  Communications 

The  nonlinear  codes  developed  for  use  in  nano-structured  storage  application  proved  to  be  valuable  for 
mitigating  the  distortion  associated  with  nonlinear  fiber  propagation.  We  have  considered  a  single 
wavelength  channel  in  a  nonlinear  dispersive  single-mode  fiber  in  the  normal  dispersion  regime.  We  study 
the  propagation  of  binary  sequences  using  OOK  RZ  modulation  at  a  bit  rate  of  10  Gb/s,  under  different 
dispersion  and  input  power  conditions.  We  determine  the  impact  that  dispersion  and  nonlinearity  have  on 
the  pulse  sequences  by  measuring  the  Euclidian  distances  among  all  pairs  of  received  sequences.  We  then 
model  the  set  of  Euclidian  distances  as  a  fully  connected  graph  and  by  means  of  a  clique-finding  algorithm 
we  search  for  subsets  of  sequences  that  form  codes  matched  to  this  nonlinear  channel,  resulting  in 
significant  bit  error  rate  improvement.  We  present  a  numerical  analysis  of  the  effects  that  GVD  and  SPM 
have  on  the  Euclidian  distances  among  pulse  sequences  at  the  receiver  in  a  single  channel  of  a  WDM 
system.  We  observe  that  nonlinearity  induces  a  decrease  in  the  minimum  Euclidian  distance  within  the 
normal  dispersion  regime  and  that  this  decrease  is  stronger  at  high  power  levels.  We  model  the  set  of 
Euclidian  distances  as  a  fully  connected  graph,  where  each  sequence  is  a  node  and  the  Euclidian  distance 
between  each  pair  of  sequences  labels  the  corresponding  edge.  We  then  remove  all  edges  with  a  Euclidian 
distance  smaller  than  an  arbitrarily  chosen  value.  Using  a  Clique  finding  algorithm,  we  search  for  subsets 
of  nodes  that  form  fully  connected  subgraphs  of  this  pruned  graph.  From  these  subsets  we  obtain  nonlinear 
codes  that  are  well  matched  to  this  nonlinear  fiber  channel,  giving  a  significant  reduction  in  bit  error  rate 
compared  to  linear  codes  of  the  same  length.  More  details  about  this  work  can  be  found  in  Appendix  2. 

V.  Information  Propagation  in  a  Fast-Light  Medium 

Projects  II  and  III  are  concerned  with  the  use  of  nanostructures  to  control  light  in  the  space-domain.  It  is 
also  possible  to  employ  nanostructures  to  control  light  in  the  time-domain.  One  manifestation  of  this 
control  is  the  engineering  of  materials  with  large  anomalous  dispersion.  The  special  theory  of  relativity 
states  that  the  velocity  of  information  propagation  is  limited  to  the  speed  of  light  in  vacuum  (c).  It  is  known 


however,  that  optical  pulses  traveling  in  a  medium  with  large  anomalous  dispersion  have  faster-than-e 
group  velocities.  The  key  to  resolving  this  controversy  is  to  quantify  how  information  is  (a)  encoded  on  an 
optical  pulse  and  (b)  affected  by  the  propagation  medium.  We  have  transmitted  two  distinct  symbols 
through  a  ’fast  light'  medium  and  quantify  both  the  propagation  time  and  the  detection  latency.  The  former 
allows  us  to  bound  the  information  velocity;  while,  the  latter  reconciles  the  apparent  discrepancy  with 
group  velocity.  In  our  experiment  we  find  that  although  the  group  velocity  vastly  exceeds  c\  the  time 
required  to  detect  information  propagating  through  the  medium  is  slightly  longer  than  the  time  required  to 
detect  the  same  information  travelling  through  vacuum.  Our  work  highlights  the  issues  that  must  be 
addressed  in  designing  devices  that  operate  at  the  ultimate  limit  imposed  by  the  special  theory  of  relativity. 
More  details  about  this  work  can  be  found  in  Appendix  3  and  related  work  on  causality  in  metamaterials  in 
Appendix  6. 

VI.  Information  Propagation  in  a  Slow-Light  Medium 

Through  manipulation/engineering  of  materials  at  the  nano-scale  it  is  possible  to  create  an  optical  medium 
in  which  the  group  velocity  (vg)  is  a  small  fraction  of  c.  These  systems  represent  a  type  of  short-term 
memory  that  can  be  valuable  in  optical  communication,  computing,  and  processing  applications.  In  our 
experiments  we  have  vy  =  c/100.  It  is  often  assumed  that  the  velocity  of  information  propagation  in  such  a 
medium  is  given  by  the  group  velocity  whenever  vg<c.  During  this  part  of  our  year  3  activities  we  sought  to 
experimentally  measure  the  velocity  of  information  propagation  in  such  a  medium.  By  imposing  a  point  of 
non-analyticity  on  the  propagating  waveform,  we  create  a  well-defined  time  at  which  information  is 
launched.  By  estimating  the  time  at  which  this  information  can  first  be  measured  we  are  able  to  obtain  an 
estimate  of  the  information  velocity  (V;).  Figure  1  shows  some  experimental  results.  Visually  we  can  see 
that  the  time  at  which  information  (i.e.,  the  discontinuity)  appears  at  the  detector  is  nearly  the  same  for  the 
signal  propagating  through  vacuum  and  the  slow-light  medium.  A  bit  error  rate  analysis  of  this  data  reveals 
that  v,  is  approximately  0.6c.  This  surprising  result  reveals  that  (a)  the  group  velocity  is  not  always  a  good 
measure  of  information  velocity,  (b)  information  can  be  carried  in  the  pre-cursor  fields  and  propagates  at 
nearly  c,  and  (c)  these  pre-cursor  fields  can  be  large  enough  to  provide  reliable  measurements  of  the 
associated  information.  We  are  still  working  to  understand  the  implications  of  these  results;  however,  it  is 
already  clear  that  the  design  of  optical  pulse  shapes  should  be  conducted  with  awareness  of  the  nature  of 
information  propagation  in  the  system  under  study.  More  details  about  this  work  can  be  found  in  Appendix 
4. 

VII.  Distortion  Compensation  in  a  Slow-Light  Delay  Line 

A  significant  limitation  to  current  methods  of  achieving  slow-light  pulse  delay  arises  from  the  distortion 
that  accompanies  this  delay.  We  have  developed  a  new  method  through  which  large  pulse  delay  may  be 
achieved  under  the  constraint  that  the  distortion  does  not  exceed  a  particular  limit.  Our  technique  is  based 
on  combining  multiple  dispersive  devices  in  order  to  engineer  the  dispersion  of  the  overall  system.  We 
have  analyzed  and  experimentally  demonstrated  the  case  of  a  gain-doublet  in  which  a  pair  of  Lorentzian 
gain  lines  are  used  to  achieve  significant  pulse  delay  relative  to  a  single-line  system.  A  factor  of  6.25 
improvement  in  delay  and  a  factor  of  2  improvement  in  pulse  bandwidth  can  be  achieved  by  use  of  our  new 
method  of  jointly  optimizing  the  overall  system  parameters.  The  optimization  of  more  complex  distortion- 
compensated  systems  based  on  resonators,  interferometers,  and  PBGs  is  currently  underway.  More  details 
about  this  work  can  be  found  in  Appendix  5. 

VIIT.  Pulse  Shaping  from  a  Spatially  Varying  Group  Index 

The  new  technique  that  we’ve  developed  is  a  spatio-temporal  technique.  It  involves  the  spatial  modulation 
of  the  pulse  delay  created  by  an  optical  medium,  as  shown  in  Figure  I .  As  the  pulse  propagates  through  the 
medium,  some  spatial  components  of  the  pulse  are  delayed  more  than  others.  This  pulse  “shear”  can  be 
used  to  control  the  temporal  pulse  shape.  As  the  pulse  is  recombined^  for  example  by  refocusing  or 
injecting  the  pulse  into  a  fiber,  the  temporal  intensity  profile  is  determined  by  the  spatial  pulse  delay  profile 
of  the  medium.  As  in  other  spatio-temporal  techniques,  the  pulse  delay  profile  need  only  be  controlled  on 
the  timescale  that  a  single  shape  is  used.  This  control  can  either  be  optical  or  electronic  depending  on  the 


pulse  delay  mechanism  employed.  Advantages  of  this  technique  over  other  spatio-temporal  techniques  are 
that  it  may  be  more  easily  extended  to  a  wide  range  of  pulse  lengths  (using  media  with  different  pulse 
delays)  and  it  can  be  very  compact.  Figure  2  shows  preliminary'  results  demonstrating  simple  pulse-shaping 
of  a  pulse  propagating  through  a  BaTi03  crystal  with  two  regions  of  different  pulse  delays.  These  regions 
were  created  by  controlling  the  intensity  of  a  pump  beam  used  to  generate  the  pulse  delay  in  each  region. 
Because  the  response  time  of  BaTi03  is  very  slow,  this  material  is  best  suited  for  shaping  long  pulses, 
whereas  short  pulses  can  be  shaped  using  the  same  technique  with  faster-response-time  materials  such  as 
nonlinear  crystals,  atomic  systems,  or  “dispersion  engineered”  materials,  possibly  allowing  electronic 
control  with  the  use  of  EO  inclusions. 

The  use  of  a  medium  with  spatially  varying  optical  properties  also  raises  the  possibility  of  diffraction.  To 
our  knowledge,  the  effect  of  diffraction  from  spatial  varying  group  delay  has  not  been  previously  studied. 
We  are  continuing  both  theoretical  and  experimental  activities  on  this  topic. 
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Figure  1 :  Medium  with  a  spatially  varying  group  index  can  be  used  as  a  pulse  shaper. 
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Figure  2:  Initial  results  demonstrating  the  operation  of  our  pulse  shaper. 


IX.  Gaussian  Beam  Scattering  from  Nanostructures 


We  have  continued  to  examine  the  scattering  from  nanostructures  in  the  presence  of  complex  material 
environments.  Particular  emphasis  has  been  given  to  nanorods  and  nanospheres  with  plasmonic  coatings. 
The  enhanced  scattering  from  these  plasmonic  nanostructures  provides  enhanced  discrimination  between 
scattering  states  for  our  optical  data  storage  application.  This  effect  has  been  used  effectively  for  our  nano¬ 
voxel  studies.  Additionally,  it  also  allows  one  to  guide  light  along  subwavelength  dimensions  to  produce 
subwavelength  field  localization.  An  example  is  shown  below.  A  transmission  line  of  nano-rods  provides 
a  means  of  capturing  light  and  guiding  it.  The  output  spot  is  significantly  subwavelength.  More  details 
about  this  work  can  be  found  in  Appendix  7.  This  was  a  paper  presented  at  a  Special  Session  on  Advances 
in  Plasmonic  and  Other  Electromagnetic  Phenomena  Co-Organized  by  R.  W.  Ziolkowski  and  N.  Engheta 
and  that  appears  in  the  Proceedings  of  the  IEEE  International  Symposium  on  Antennas  and  Propagation, 
Washington  DC,  July  2005 


Figure  3:  A  nano-rod  array  transmission  line  is  used  to  guide  and  to  localize  500nm  light  acquired  from  an 
incident  Gaussian  beam  to  sub-wavelength  resolutions. 
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Abstract 

A  new  optical  data  storage  structure  is  proposed  that  utilizes  surface  relief  profiles  at  sub-wavelength  scales  to  dis¬ 
tinguish  memory  slates  in  a  read-only  memory  system.  Groupings  of  power  sensors  directly  beneath  the  features  read 
light  transmitted  through  the  structure  from  above.  These  are  used  to  distinguish  states  from  one  another,  where  dif¬ 
ferent  surface  profiles  lead  to  different  states.  Two-dimensional  simulations  are  performed  using  the  finite-dilTerence 
time-domain  (FDTD)  method  to  optimize  parameters  of  the  system,  and  to  determine  the  data  capacity  and  density 
of  optimal  systems.  Two  optimal  systems  are  studied  in  detail;  one  intended  as  a  polymer  construction,  the  other  as 
silicon.  Techniques  are  developed  to  estimate  capacity  and  density  from  subsets  of  simulations.  Simulations  estimate 
that  densities  of  3.17  GBits/crrr  are  possible  with  this  type  of  system. 

©  2005  Elsevier  B.V.  All  lights  reserved. 

PACS:  42.79. Vb 

Keywords:  Optical  storage;  Memory  systems 


1.  Motivation 

Any  system  with  a  large  number  of  distinguish¬ 
able  states  can  act  as  a  memory  (storage)  device. 
Current  generation  high-capacity  memory  devices 


Corresponding  author.  Tel.;  +1  403  220  4129;  fax:  +1  403 
282  6855. 

E-mail  address:  mpotter@ucaIgary.ca  (M.E.  Potter). 


(archival  rather  than  e.g.  flash)  generally  rely  on 
one  of  two  technologies:  optical  or  magnetic.  Opti¬ 
cal  systems  include  CD-ROMs,  DVDs,  and  the 
next  generation  of  Blu-Ray  discs.  For  all  three, 
states  are  created  and  distinguished  from  each 
other  by  surface  features  (pits  and  ridges);  each 
feature  represents  a  single  binary  state.  CDs  are 
limited  by  classical  resolution  limits  of  approxi¬ 
mately  one  bit  per  square  wavelength.  DVDs  and 
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Blu-Ray  discs  use  multiple  layers  and  can  achieve 
four  to  eight  bits  per  square  wavelength.  In  the 
case  of  magnetic  storage  media,  data  is  encoded 
in  the  magnetic  polarity  of  small  sub-domains  of 
the  medium.  Scanning  a  probe  over  the  structure 
reads  the  polarity  in  given  regions.  Demonstra¬ 
tions  by  Kong  et  al.  [1]  have  shown  magnetic  stor¬ 
age  densities  on  the  order  of  10  Gbits/in2 
( 1 .55  Gbits/cm2).  Sub-wavelength  topographical 
features  created  with  lithography  can  also  be  used 
to  imprint  data  [2].  Features  are  then  probed  using 
atomic  force  microscopy  techniques,  since  lasers 
cannot  be  effectively  used  to  probe  each  individual 
feature.  Impressive  storage  densities  on  the  order 
of  400  Gbits/in2  (62  GBits/cm2)  have  been  demon¬ 
strated  [3].  One  major  limitation  with  disc-based 
systems  is  that  each  feature  is  read  serially,  and 
the  time-scale  of  the  reading  process  introduces  la¬ 
tency  effects. 

An  alternate  optical  storage  methodology  that 
shows  promise  is  holography.  The  advantage  of 
holography  is  that  multiple  data  “pages5'  can  be 
stored  in  the  same  medium  simply  by  altering  the 
angle  at  which  the  hologram  is  recorded  in  the 
medium.  This  is  a  form  of  parallel  data  addressing 
since  an  entire  page  can  be  recovered  at  once. 
Impressive  storage  densities  of  10-40  Gbits/cm2 
have  been  predicted  for  such  systems,  and  demon¬ 
strations  with  holographic  RAM  have  indicated 
faster  access  times  than  hard  discs  for  computers 

The  focus  of  our  research  is  to  examine  systems 
that  will  capitalize  on  parallel,  near-field  readout. 
Our  intended  application  for  this  system  would 
be  analogous  with  CD-ROM  or  DVD  uses,  or 
any  usage  that  requires  multiple  replicas  of  the 
same  data.  A  target  application  would  be  for  sys¬ 
tems  requiring  low  latency,  such  as  video  games. 
Other  applications  would  be  for  mass  distribution 
of  high  content  media  such  as  movies,  games,  and 
reference  materials.  Specifically,  we  postulate  a 
read-only  memory  system  that  rivals  magnetic 
storage  technologies  for  capacity/density,  but  uti¬ 
lizing  optical  readout  at  a  wavelength  of 
A  =  500  nm.  States  will  be  defined  by  groupings 
of  sub-wavelength  surface  features,  which  are 
identified  by  reading  the  signal  (in  the  near-field) 
transmitted  through  the  structure.  The  aim  is  to 


achieve  a  density  of  16  Gbits/cm2  in  order  to  ap¬ 
proach  the  same  densities  as  magnetic,  while  sur¬ 
passing  DVD  and  Blu-Ray  densities.  This  paper 
outlines  the  details  of  the  system  in  Section  2, 
and  the  modeling  technique  used  is  introduced  in 
Section  3.  Results  from  initial  simulations,  discus¬ 
sion,  and  conclusions  are  then  presented  in  Sec¬ 
tions  4~6. 


2.  System  description 

2. 1 .  Intended  implementation 

The  proposed  system  consists  of  a  series  of 
wavelength-scale  nanostructures  (cells),  with  de¬ 
signed  sub-wavelength  features  which  will  be  dis¬ 
tinguished  using  near-field  optical  measurements. 
The  features  are  formed  by  introducing  dielectric 
inhomogeneity.  For  the  purposes  of  this  study, 
the  system  is  considered  in  two  dimensions.  The 
entire  system  consists  of  a  row  of  cells  on  a  flat  sur¬ 
face,  of  which  three  cells  are  shown  in  Fig.  1  (in  the 
third  dimension,  multiple  rows  would  be  added). 
The  entire  system  is  illuminated  from  above  by 
an  optical  source  at  a  wavelength  of  A  =  500  nm. 
Sensors  directly  underneath  each  cell  record  opti¬ 
cal  power  measurements  in  parallel.  These  mea¬ 
surements  are  used  to  identify  the  profile  of  the 
sub-wavelength  features,  which  is  then  mapped 
to  the  stream  of  bits  it  is  meant  to  represent.  Since 
the  entire  system  is  illuminated  at  once,  the  system 
does  not  require  moving  components  and  the  data 
stored  in  the  system  is  read  in  parallel  at  the  phys¬ 
ical  level. 

Each  cell  is  one  wavelength  wide,  although  the 
entire  system  will  consist  of  many  such  cells  on 
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Fig.  I.  System  layout. 


58 


M.E.  Potter  ct  ai  I  Optics  Communications  253  (2005  )  56-69 


something  similar  to  a  compact  disc.  For  the  pur¬ 
poses  of  analysis,  we  focus  on  one  cell  of  interest 
(COI)  and  two  neighbouring  cells.  This  allows  us 
to  analyze  the  uncertainty  of  the  measurements 
introduced  by  the  proximity  of  neighbouring  cells 
via  crosstalk  phenomena.  All  the  cells  in  a  system 
have  the  same  dimensions  as  described  by  the  set 
of  parameters  indicated  in  Fig.  2.  Three  horizontal 
layers,  namely  the  fill,  feature,  and  base  layers, 
make  up  a  cell.  The  fill  and  base  layers  are  homo¬ 
geneous  materials.  The  feature  layer  is  divided  into 
different  sections  allowing  for  either  the  fill  or  base 
material  to  be  present,  creating  a  non-homoge- 
neous  layer  and  therefore  the  desired  sub-wave¬ 
length  features. 

The  system  is  intended  to  be  used  as  a  ROM  de¬ 
vice,  whose  low  latency  and  high  density  may  jus¬ 
tify  the  cost  premium  associated  with  silicon 
fabrication  processes.  In  fact,  we  anticipate  that 
the  surface  relief  structures  can  be  stamped,  and 
thus  may  benefit  from  the  economies  of  scale  pos¬ 
sible  from  mass  production  technologies  such  as 
for  CD-ROMs  and  DVDs. 

We  consider  that  the  fabrication  of  the  detector 
layer  wil  l  be  asynchronous  with  that  of  the  feature 
layers,  and  as  such  there  will  have  to  be  a  calibra¬ 
tion  procedure  to  tune  for  alignment  tolerance. 
For  instance,  a  set  of  calibration  elements  would 
be  fabricated  at  higher  tolerances  and  analyzed. 
Each  storage  “disc”  would  contain  a  few  fiducial 
elements,  and  the  readings  from  these  compared 
to  the  calibrated  standards  would  allow  us  to  select 
the  correct  decoder  to  extract  the  stored  data  on 
the  rest  of  the  “disc.”  The  detectors  themselves  will 
not  require  individual  sense  amps,  but  rather  we 
anticipate  a  design  that  would  be  a  matrix- 
addressed  photo-conductor  array  requiring  signal 
conditioning  only  for  each  row  and  column. 


2.2.  Distinguishing  memory  states  and  capacity 

The  current  research  focuses  on  determining  the 
capacity  for  a  two-dimensional  nanostructured 
system  with  the  aim  to  map  as  many  features  per 
COI  as  possible  (/if).  Thus,  for  example,  with 
Hf  =  5  features  per  cell  there  are  23  =  32  different 
possible  COI  configurations  (states).  Parameteri¬ 
zation  was  performed  to  find  systems  that  allow 
maximal  distinguishability  amongst  the  states. 
Parameters  that  were  varied  include;  the  feature 
height,  the  base  height,  the  gap  size,  and  the  index 
of  refraction  of  the  base  and  feature  layers. 

Different  states  are  distinguished  using  near-field 
optical  power  measurements  from  sensors  located 
directly  underneath  each  COI.  Each  sensor  pro¬ 
duces  a  power  measurement,  P,,  and  each  cell  has 
an  array  of power  sensors  resulting  in  a  power  vec¬ 
tor,  P  =  (Pi ,  Po,  • .  • .  P*f).  By  altering  the  configura¬ 
tion  of  features  for  a  given  COI,  there  are  2',f 
possible  power  vectors  that  may  be  used  to  represent 
different  states.  Note  that  we  do  not  seek  a  one-to- 
one  mapping  between  features  and  representative 
bits:  any  state  vector  may  be  used  to  represent  any 
bit  stream  we  wish.  In  spite  of  the  desired  scatter 
from  features  within  a  cell,  there  will  also  be  undesir¬ 
able  scattering  from  adjacent  cell  structures  (neigh¬ 
bours)  which  will  introduce  crosstalk.  This  crosstalk 
will  manifest  itself  in  uncertainty  in  the  measure¬ 
ment  of  the  state  vector  for  a  given  cell.  This  uncer¬ 
tainty  may  mean  that  some  state  vectors  could 
become  indistinguishable  and,  hence,  they  may 
not  be  used.  Thus,  the  loss  of  distinguishability  re¬ 
duces  the  capacity  of  that  system  parameterization. 
Obviously  the  goal  is  to  find  the  system  parameter¬ 
ization  that  provides  the  maximum  capacity.  This  is 
accomplished  by  finding  all  the  state  vectors  that  are 
distinguishable  from  all  other  state  vectors. 
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Fig.  2.  Geometric  parameters  of  the  system. 
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It  should  be  repeated  that  the  capacity  is  limited 
only  by  the  number  of  distinguishable  patterns 
that  can  exist  in  the  multi-dimensional  space  of 
the  power  vector.  This  allows  us  to  take  advantage 
of  many  other  degrees  of  freedom  such  as  source 
wavelength,  feature  heights,  etc.  For  example,  the 
height  of  the  feature  above  each  detector  may  be 
four-tiered  instead  of  two-tiered,  giving  the  possi¬ 
bility  of  4 5  =  1024  possible  states,  and  hence  the 
potential  for  more  than  one  bit  per  detector. 

3.  Modeling  methodology 

3.1.  Methodology 

To  investigate  different  parametenzations  of  our 
proposed  system,  it  was  necessary  to  simulate  differ¬ 
ent  systems  with  changing  physical  dimensions  and 
dielectric  properties.  To  meet  our  goal  of  16  Gbits/ 
cm2,  a  system  with  5-7  features/cell  is  required. 
Even  at  the  lower  limit  of  5  features/cell,  to  obtain 
an  exhaustive  characterization  (to  examine  cross¬ 
talk  for  a  single  system)  would  require  us  to  simulate 
each  possible  COI  configuration  (25  =  32),  and  for 
every  possible  neighbour  configuration  (2 10  = 
1024)  requires  approximately  32,000  simulations. 
Fully  investigating  several  different  systems  quickly 
becomes  prohibitively  time-consuming. 

In  the  first  stage  of  our  methodology,  COIs  are 
analyzed  and  parameterized  in  isolation.  Then  full 
sets  of  simulation  for  systems  with  3  features/cell 
(512  simulations/system)  were  performed  to  see  if 
we  could  extrapolate  what  fraction  of  the  five- 
feature  simulations  would  be  necessary  for  reason¬ 
able  estimates  of  capacity.  This  analysis  method  is 
described  in  further  detail  in  Section  4.  A  further 
benefit  of  this  three-feature  analysis  is  that  it  al¬ 
lows  us  to  visualize  the  clustering  of  power  mea¬ 
surements  in  three-dimensional  visual  mappings. 

3.2.  Modeling  tool 

The  modeling  tool  must  be  well-developed  and 
able  to  handle  complex  geometries,  generating 
reliable  results  so  that  we  can  concentrate  on  the 
analysis  and  statistics.  For  these  reasons  a  two- 
dimensional  (TM  polarization)  finite-difference 


time-domain  (FDTD)  method  was  chosen  due  to 
its  flexibility  in  materials  and  geometries.  It  is  a 
reliable  and  time-proven  method,  and  is  particu¬ 
larly  well  suited  for  heterogeneous  structures. 
Maxwell’s  equations  in  differential  form  are  dis¬ 
cretized  in  time  and  space,  and  explicit  expressions 
determine  the  electric  and  magnetic  fields  at  future 
time  steps  based  solely  on  their  values  at  previous 
time  steps  and  on  the  local  electromagnetic  proper¬ 
ties  of  materials.  More  details  can  be  found  in  [7]. 

The  simulation  space  was  discretized  at  a  reso¬ 
lution  of  A*  =  Az  =  5  nm  (100  cells  per  free  space 
wavelength)  which  because  of  stability  condi¬ 
tions  results  in  a  time  discretization  of 
At  -  1 1.55  as.  The  absorbing  boundary  condition 
at  the  edge  of  the  simulation  space  was  a 
L2TDLM  (see  [8])  of  1 1  layers,  and  the  simula¬ 
tions  were  allowed  to  progress  for  3000  timesteps 
-  approximately  20  cycles  of  the  source. 

The  optical  source  was  a  continuous  wave 
Gaussian  beam  with  TM  polarization  and  a  wave¬ 
length  of  a  =  500  nm.  For  simulation  purposes,  the 
beam  is  generated  by  invoking  the  equivalence 
principle  at  a  line  lying  along  the  waist  of  the  beam 
(see,  for  example,  the  total/scattered  field  formula¬ 
tion  in  [7]).  The  waist  of  the  beam  was  three  wave¬ 
lengths  in  width  in  our  simulations,  and  was 
approximately  two  wavelengths  above  the  struc¬ 
ture.  Observation  tasks  consisted  of  producing 
time  histories  of  the  electric  and  magnetic  fields 
at  each  of  the  sensor  locations,  from  which  the 
Poynting  vector  and  the  associated  power  mea¬ 
surements  were  obtained  in  a  post-processing 
phase.  Note  that  reported  power  measures  are  in 
arbitrary  units,  which  will  scale  accordingly  as 
the  fields  are  scaled.  The  source  used  is  identical 
for  all  simulations  reported,  so  that  comparisons 
can  be  made.  In  this  sense,  wre  are  simulating  that 
one  detector  lies  directly  beneath  each  feature  in 
the  COI  (i.e.  one  detector  per  feature),  and  that 
the  fill-factor  of  the  detector  is  100%  that  is, 
we  assume  that  each  detector  completely  covers 
the  plane  beneath  each  feature. 

3.3.  Structures  modeled 

Two  possibilities  for  systems  were  envisaged:  one 
based  on  polymers  that  can  be  easily  stamped  to 
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create  the  features,  and  one  based  on  silicon  for 
which  lithographic  processes  will  probably  be  neces¬ 
sary  to  create  the  features.  To  account  for  these  two 
cases  in  the  FDTD  simulations,  only  two  refractive 
index  profiles  were  used.  For  the  former,  the  base 
and  feature  layers  had  an  index  of  refraction  of 
polymer  =  1  -8;  and  for  the  latter  the  base  and  feature 
layers  had  an  index  of  refraction  of  /^silicon  -  3.5.  In 
both  cases  a  fill  layer  with  an  index  of  refraction  of 
«r,i]  =  1 .2  was  used,  as  well  as  a  substrate  with  an  in¬ 
dex  of  refraction  of  flsub»trat«  =  3-5. 


4.  Analysis  and  results 

To  analyze  capacity,  it  is  necessary  to  determine 
the  number  of  distinguishable  states  in  a  given  sys¬ 
tem.  For  a  given  COI,  reflecting  and  diffracting 
fields  from  neighbouring  cells  interact  to  produce 
interference  in  the  COI  power  measurement.  If 
neighbouring  cells  were  infinitely  far  away,  then 
there  would  be  no  interference  and  no  ambiguity 
in  the  assignment  of  a  particular  state  to  a  partic¬ 
ular  power  measure  (Fig.  3(a)).  When  neighbour¬ 
ing  cells  are  close,  the  effect  of  the  interference  is 
to  cause  a  slight  shift  in  the  measured  power  vector 
(Fig.  3(b)),  and  uncertainty  is  introduced  in  assign¬ 
ment  of  the  state  since  different  neighbour  config¬ 
urations  will  affect  it  to  differing  degrees.  When 
the  regions  of  uncertainties  for  two  different  states 
are  close,  there  is  a  greater  chance  for  mis-assign- 
ment  of  a  state.  The  goal  of  our  analysis  was  to 
establish  how  many  states  for  a  given  system  are 


distinguishable,  and  hence  determine  the  system 
capacity. 

Analysis  consisted  of  an  initial  analysis  to  deter¬ 
mine  systems  that  provide  maximum  separation  of 
state  vectors  in  isolation  (i.e.  no  crosstalk).  Good 
candidates  were  then  rigorously  analyzed  with 
neighbouring  cells,  to  determine  measures  of 
uncertainty  and  thus  capacity. 

4.1.  Initial  analysis 


In  the  first  stage,  COls  were  analyzed  in  isolation 
to  get  a  initial  estimation  of  which  parameters  pro¬ 
vide  the  greatest  distinguishability  in  state  vectors. 
Because  there  are  no  neighbours,  all  state  vectors 
are  unique.  Consequently,  the  aim  was  to  identify 
parameterizations  that  provide  maximum  separa¬ 
tion  (Euclidean)  distance  between  state  vectors. 
Intuitively  we  should  expect  that  systems  with  more 
separation  will  provide  better  distinguishability. 

To  understand  the  metrics  we  define  for  evaluat¬ 
ing  systems,  first  consider  a  set  of  n  power  vectors, 
with  an  individual  power  vector  being  described  by 

P  i  =  {Pn,Pa>--,Pin(l  1 

Then  we  define  the  vector  describing  the  cen¬ 
troid  of  the  distribution 
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The  distance  from  one  vector  to  another  vector  is 
(in  the  Euclidean  sense)  dfj  =  IP,*  -  P;|,  and  the  dis¬ 
tance  from  a  vector  to  the  centroid  is  Jic  =  (P,-  —  Pc|. 


f  P2 


State  1 


Fig.  3,  The  effect  of  interference  from  neighbouring  COIs  for  a  simplified  system  (two  features).  In  (a),  for  COls  in  isolation  (i.e. 
neighbouring  cells  at  infinity)  would  exhibit  no  uncertainty  in  identifying  two  different  stales  power  measure  would  be  exact.  In  (b), 
neighbouring  cells  are  close,  and  for  a  given  COI  state,  changes  in  the  neighbouring  cells  will  cause  interference  and  hence  uncertainly 
in  the  slate  power  measurement. 
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The  average  distance  from  the  centre  is  then  given 
by  dc  c-  We  define  several  variables  to  be 

able  to  assess  the  merit  of  a  system: 

diststd  —  ^  y/^l^  —  dc)1  the  standard 
deviation  in  distance  of  vectors  from  the  aver¬ 
age  vector  location, 

avgpdist  -  { d,j )  -  (|P,-  —  Pj\)  the  average  in  dis¬ 
tance  between  vectors, _ 

pdiststd  =  =  yj( |  P,  -  P7|2)  the  standard 

deviation  in  distance  between  vectors. 

The  average  and  standard  deviation  in  power 
for  each  detector  were  also  identified  and  compared 
graphically.  Again,  if  the  power  measurements 
were  distributed  more  widely  throughout  five- 
space,  we  would  expect  better  distinguishability. 

For  the  geometry  of  the  structures  in  the  simu¬ 
lations,  the  fill  depth  was  kept  constant  at 
f/ini  =  100  nm,  and  the  width  of  each  feature  was 
always  Hycaturc  ~  100  nm.  For  parameterization, 
the  height  of  the  features  is  varied  between 
h  =  50  and  h  =  250  nm,  and  the  depth  of  the  base 
layer  is  varied  between  d  -  10  and  d  -  500  nm. 

Results  are  found  in  Table  1.  In  general,  the 
higher  refractive  index  contrast  provides  better 
separation  between  state  vectors.  This  can  be  seen 
better  by  observing  Fig.  4.  The  power  at  each 
detector  for  all  32  configurations  is  plotted  as 
points,  the  mean  power  at  each  detector  is  plotted 
as  the  solid  line,  and  error  bars  indicate  two  stan¬ 
dard  deviations  in  the  measured  power  from  the 
mean  at  each  detector.  Examples  of  poor  systems 
and  good  systems  are  shown  for  both  cases  of 
refractive  index  profiles.  Though  there  are 
probably  several  choices  of  systems  that  could  be 
made  for  more  detailed  studies,  we  will  focus  on 
two  -  one  n  -  1.8  (polymer)  and  one  n  ~  3.5  (sili¬ 
con)  -  that  show  the  best  promise.  Systems  are  des¬ 
ignated  with  the  refractive  index  /?,  base  depth  d , 
and  feature  height  /?,  as  nXXdYYhZZ.  Using  this 
naming  convention,  the  chosen  best  polymer  sys¬ 
tem  is  /7l.8r/500/7250  since  it  excels  in  all  three  mea¬ 
sures  compared  to  the  others.  The  chosen  best 
silicon  system  is  //3.5c/l 0/z  1 50;  although  the  n3. 
5d5Q0h] 50  is  similar  in  measures,  we  have  opted 
to  look  at  the  former  because  it  would  require  less 


material  for  fabrication.  The  chosen  systems  are 
indicated  in  Fig.  5,  and  are  the  basis  for  all  re¬ 
ported  results  from  this  point  on. 

4.2.  Distinguishable  states  and  capacity  measures 

To  determine  how  well  candidate  systems  per¬ 
form  in  the  presence  of  crosstalk  (i.e.  with  neigh¬ 
bouring  cells),  we  first  had  to  define  metrics  for 
this  crosstalk.  With  a  full  set  of  simulations  for 
every  possible  COI  and  neighbour  configuration, 
this  allows  us  to  identify  distinguishable  states, 
and  therefore  capacity. 

For  the  purposes  of  this  study,  we  designate  a 
region  of  uncertainty  for  a  given  COI  by  a  power 
vector  cloud,  where  all  of  the  possible  interference 
conditions  lie.  This  region  is  then  designated  by  a 
centre  C  and  radius  Ry  as  demonstrated  in  Fig,  6 
for  an  example  two-feature  system.  The  centre  is 
defined  by  the  arithmetic  mean  of  all  the  possible 
vectors  in  the  region,  and  the  radius  is  defined  as 
the  scalar  distance  between  this  centre  and  the  fur¬ 
thest  outlying  measurement.  For  the  three-feature 
system,  regions  are  described  by  spheres,  and  be¬ 
yond  three  features  by  hyperspheres  of  the  appro¬ 
priate  dimension. 

For  a  COI  to  be  distinguishable,  it  is  necessary 
that  its  region  of  uncertainty  be  distinct  from  all 
other  regions  for  all  other  COls,  otherwise  those 
states  can  not  be  used.  This  allows  us  to  determine 
the  capacity  of  the  system.  As  an  example,  in  Fig.  7 
a  two-feature  system  is  shown  with  the  regions  of 
uncertainty  for  each  of  the  four  possible  COI  con¬ 
figurations.  In  this  case,  two  of  the  states  are  indis¬ 
tinguishable  from  each  other  in  the  presence  of 
crosstalk  and  therefore  can  not  be  used.  Thus  this 
system  has  a  storage  capacity  of  2  states/COI 
(1  bit/COI). 

In  point  of  fact,  this  method  of  determining 
capacity  assumes  that  if  two  states  are  overlap¬ 
ping,  then  they  are  both  discarded.  In  reality,  we 
only  need  discard  one  of  the  two.  The  choice  of 
which  state  to  discard  is  complicated  however,  be¬ 
cause  this  depends  on  its  relationship  with  all  of 
the  other  states.  The  choice  of  which  states  to  keep 
and  which  to  discard  is  actually  an  NP-hard  prob¬ 
lem.  For  instance,  even  with  32  possible  states  for 
our  COI,  there  are  232  possible  subsets  of  states 
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Table  l 


Initial  analysis  of  several  systems  in  isolation 


Refractive  index  of  base  layer 

Depth  of  base  4  (nm) 

Height  of  feature  h  (nm) 

distsui 

avgpdist 

pdiststd 

1.8 

50 

10.85 

120.18 

30.55 

10 

150 

46.16 

286.96 

86.70 

250 

150.54 

577.22 

223.86 

50 

11.47 

80.55 

22.71 

100 

150 

120.28 

397.07 

171.30 

250 

186.69 

598.97 

264.03 

50 

4.80 

139.74 

70.47 

500 

150 

137,94 

426.24 

210.88 

250 

198.38 

777.97 

352.26 

3.5 

50 

58.71 

21 1.09 

85.31 

10 

150 

155.07 

848.95 

275.54 

250 

247.19 

583.26 

310.94 

50 

156.19 

469.51 

230.78 

100 

150 

230.63 

753.23 

330.48 

250 

153.69 

627.51 

228.06 

50 

133.01 

560.83 

233.23 

500 

150 

206.51 

836.29 

285.81 

250 

253.70 

782.23 

342.44 

that  would  need  investigating.  We  are  currently 
investigating  optimization  and  clique-finding 
methods  for  this  NP-hard  problem,  to  be  reported 
at  a  future  date. 

For  our  system  to  rival  current  magnetic  stor¬ 
age  densities,  we  require  at  least  5  features  per 
wavelength.  As  mentioned  before,  an  exhaustive 
set  of  simulations  for  one  such  system  would  re¬ 
quire  over  30,000  simulations,  which  with  the  gi¬ 
ven  choice  of  simulation  parameters  and  the 
available  computer  resources  would  take  several 
months  to  complete.  In  order  to  reduce  this 
requirement,  we  assume  that  a  random  subset  of 
neighbour  configurations  will  allow  us  to  quantify 
the  system  -  the  centres,  radii,  and  hence  capacity 
-  within  a  certain  measure  of  error.  In  order  to 
determine  what  fraction  of  the  full  set  is  necessary, 
we  performed  simulations  on  the  candidate  sys¬ 
tems  with  three  features  instead  of  five  (512  instead 
of  over  30,000  simulations).  This  then  gives  us  a 
full  set  to  compare  measures  against  as  we  increase 
the  fraction  of  neighbours  used.  Our  intention  is 
then  to  use  this  knowledge  to  estimate  the  fraction 
necessary  for  five-feature  simulations. 


The  capacity  (and  density)  will  be  affected  by 
the  spacing  between  cells  (the  gap  size);  as  the 
spacing  gets  larger  and  larger,  we  should  expect 
better  and  better  distinguishability,  and  hence  bet¬ 
ter  capacity.  However,  the  density  (capacity  per 
linear  measure)  should  peak  at  some  optimal  spac¬ 
ing.  We  wished  to  make  this  spacing  a  constant  in 
all  future  simulations,  and  so  before  proceeding 
with  full  simulations  for  capacity,  we  did  some 
nominal  simulations  with  a  very  small  subset  of 
neighbours,  but  increasing  gap  size  to  get  the  den¬ 
sity  trend.  Fig.  8  shows  the  results  of  these  initial 
simulations  as  the  gap  size  increases  from  50  to 
500  nm.  We  see  that  the  turnover  point  for  capac¬ 
ity  vs.  density  lies  in  approximately  the  250- 
400  nm,  depending  on  whether  we  are  using  the 
silicon  or  polymer  system.  In  order  to  aim  for 
maximum  density,  we  thus  choose  the  250  nm 
measure  as  our  gap  size  for  all  future  simulations. 

4.3.  Three-feature  measures 

Full  sets  of  simulations  were  performed  for  the 
two  candidate  systems  indicated  in  Fig.  5,  which 
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Fig.  4.  Example  candidate  systems  in  polymer  ((a)  and  (b))  and  silicon  ((c)  and  (d)).  The  power  at  each  detector  for  all  32 
configurations  is  plotted  as  points;  mean  power  at  each  detector  is  plotted  as  a  solid  line;  error  bars  indicate  two  standard  deviations  in 
the  measured  power  from  the  mean.  Clearly  there  is  a  marked  difference  in  the  separation  of  power  measures  for  poor  and  good 
candidates. 
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Fig.  5.  Candidate  system  parameters  in  both  silicon  (left)  and  polymer  (right)  that  show  promise  for  good  capacity. 


included  the  COI  along  with  two  neighbouring 
cells.  The  cells  were  separated  from  each  other 
by  250  nm.  Other  modeling  parameters  were  as 
discussed  in  Section  3.  The  main  output  of  these 
full  sets  was  a  complete  set  of  power  vectors  for 
each  and  every  possible  configuration  of  COI 
and  neighbour. 


The  full  set  of  power  vectors  determined  our 
comparison  benchmark  for  the  actual  centres  and 
radii  of  the  power  vector  clouds.  We  then  per¬ 
formed  statistical  sub-samplings  at  certain  frac¬ 
tions  of  the  full  set  of  neighbours  (in  this  case,  64 
is  a  full  set),  and  determined  the  centre  and  radii 
with  these  random  sub-samples.  This  was 
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Fig.  6.  For  a  given  COl,  as  the  neighbour  configuration 
changes,  the  measured  power  vector  P  will  vary  slightly.  The 
variance  is  represented  by  an  uncertainty  cloud  with  a  centre  C 
and  a  radius  R. 


Fig.  7.  A  two-element  system,  indicating  a  graph  of  the  power 
vectors.  In  this  case,  COls  #2  and  #4  overlap,  and  so  are 
indistinguishable  from  each  other. 


performed  100  times  at  each  fraction  ranging  from 
3-of-64  to  64-of-64,  and  at  each  fraction  the  RMS 
error  of  the  100  random  sub-samples  (relative  to 
the  actual  radius  from  the  full  set)  was  found  for 
the  centre  and  radius.  The  results  of  this  fractional 
randomization  are  displayed  in  Figs.  9  and  10. 

Note  that,  as  expected,  as  the  fraction  taken  in¬ 
creases,  the  RMS  error  decreases  from  around 
50%  to  zero  at  the  full  fraction.  The  error  in  both 
measures  show  a  similar  trend  in  this  regard, 
although  there  is  less  error  in  the  estimate  of  the 
centre.  The  utility  in  such  a  graph  comes  from 
attempting  to  infer  what  a  necessary  fraction  is 
for  a  tolerable  level  of  error.  For  instance,  if  an  er¬ 
ror  of  10%  is  allowable  in  this  measure,  then  these 


3.5 


x  10 


Density  vs  Gap  Size 


2.5 


Q 


S  2 


ns 

U  1.5 


x 

+ 


+  + 


*  x 


+  * 
X 


1 


it  Polymer 
+  Silicon 


200  300  400  500 

Gap  Size  (nm) 


Fig.  8.  Graph  showing  the  trend  in  density  as  the  spacing 
between  cells  increases,  for  both  polymer  and  silicon  systems. 


graphs  would  indicate  that  a  fraction  of  approxi¬ 
mately  40-60%  of  the  full  set  is  required.  This 
inference  means  that,  in  the  hve-feature  case,  that 
we  should  expect  to  have  to  perform  simulations 
for  every  COI  and  a  random  subset  of  around 
400  600  of  the  1024  possible  neighbour  configura¬ 
tions.  Unfortunately  this  is  not  the  substantial 
reduction  in  the  simulation  time  that  we  were  hop¬ 
ing  for. 


RMS  Error  for  Polymer  3  Feature  System 


Fig.  9.  The  RMS  error  in  the  cenlre  and  radius  for  the  ihree- 
feature  polymer  system. 
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RMS  Error  for  Silicon  3  Feature  System 
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Fig.  10.  The  RMS  error  in  the  centre  and  radius  for  the  three- 
feature  silicon  system. 


Fraction  of  Vector  Set 

Fig.  1 1.  The  RMS  error  in  the  centre  and  radius  for  the  five- 
feature  silicon  system,  overlaid  on  the  silicon  three-feature 
system  previously  shown  in  Fig.  10. 


4.4.  Five-feature  measures 

In  order  to  verify  our  hypothesis  about  whether 
the  inference  of  the  necessary  fraction  from  the 
three-feature  studies  can  be  carried  over  to  the 
five-feature  case  (and  perhaps  to  the  //-feature 
case),  we  completed  one  exhaustive  full  set  of  sim¬ 
ulations  for  the  five-feature  silicon  case.  Simulation 
parameters  were  identical  to  the  three-feature  case, 
except  that  the  cell  of  one  wavelength  in  size  was 
divided  in  five  instead  of  three. 

Using  this  full  set  of  power  vectors,  we  then  per¬ 
formed  the  same  statistical  sub-sampling  as  previ¬ 
ously  mentioned,  over  the  range  of  fractions 
spanning  from  0  to  1.  The  results  of  this  are  shown 
in  Fig.  11,  overlaid  on  the  previous  results  for  the 
three-feature  in  silicon.  It  is  interesting  to  note  that 
in  this  case  the  error  starts  at  the  same  levels,  but 
decreases  much  more  rapidly  than  in  the  three- 
feature  case.  In  this  case,  if  an  error  of  10%  is 
allowable  we  only  need  to  take  a  fraction  of 
approximately  0.1,  corresponding  to  a  substantial 
reduction  in  simulation  time. 

Since  wc  had  a  full  set  of  data,  we  were  then 
able  to  find  our  estimate  of  the  number  of  distin¬ 
guishable  states,  and  hence  the  capacity.  For  this 
particular  system,  the  capacity  was  determined  to 
be  13  of  32  slates,  using  the  algorithm  described 
at  the  beginning  of  this  section.  The  final  check 


in  the  methodology  was  to  see  at  what  fraction 
the  capacity  is  well  predicted.  Fig.  12  shows  the 
estimated  capacity  as  the  fraction  of  neighbours 
taken  is  increased.  From  this  graph,  it  appears  that 
the  capacity  can  be  estimated  well  at  a  fraction  of 
around  0.2.  From  the  previous  graph,  this  indi¬ 
cates  that  an  error  of  around  7%  is  tolerable  for 
our  estimates  of  capacity. 

4.5.  Ellipsoid  algorithm 

One  of  the  benefits  of  performing  the  three- 
feature  simulations  is  that  we  are  able  to  visualize 
the  distribution  of  the  state  vectors  in  three-dimen¬ 
sional  space  (a  luxury  not  afforded  to  us  for  great¬ 
er  levels  of  features).  For  instance,  the  distribution 
of  the  states  (the  end-point  of  their  vectors)  for  one 
particular  system  are  indicated  in  Fig.  13.  From 
this  figure,  one  can  see  that  the  distributions  for 
particular  COIs  tend  to  be  quite  elongated  in  par¬ 
ticular  directions.  This  indicates  that  spheres  are 
not  necessarily  the  best  way  of  representing  the 
distributions,  as  they  tend  to  over-predict  the  vol¬ 
ume  occupied  by  the  distribution.  It  would  seem 
that  ellipsoids  could  be  better  predictors  of  the 
occupied  volume  of  states,  and  using  these  ellip¬ 
soids  would  provide  better  estimations  of  the 
capacities  of  systems. 
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Fig.  12.  The  estimated  capacity  as  a  function  of  the  fraction  of 
the  neighbours  taken.  At  each  fraction,  100  random  neighbours 
were  chosen  and  the  capacity  calculated.  The  error  bars  indicate 
one  standard  deviation  on  either  side  of  the  mean 


To  test  this  methodology,  we  developed  an 
algorithm  which  defines  the  enclosing  ellipsoid 
for  a  given  COPs  distribution  as  follows.  First, 
the  centre  of  the  distribution  is  found.  We  then 
find  the  furthest  outlying  state  point,  and  define 
the  vector  between  the  centre  and  this  point  as 
the  major  axis  of  the  ellipsoid.  All  state  points 
are  then  projected  onto  the  plane  perpendicular 
to  this  axis.  The  next  largest  axis  of  the  ellipse  is 
then  defined  in  a  similar  manner,  by  finding  the 
furthest  outlier  on  this  plane.  The  process  is  re¬ 
peated  until  all  axes  are  found  -  three  for  the 
three-feature  case,  or  five  for  the  five-feature  case 
(thus  defining  a  hyperellipsoid). 

From  these  defining  ellipsoids,  we  can  then 
determine  the  overlap  between  two  different  states. 
Because  the  ellipsoids  provide  a  better  indication 
of  the  states’  defining  volumes,  we  expect  the 
capacity  estimates  for  systems  to  increase  from 
previous  estimates  done  with  defining  spheres. 
For  example,  in  Figs.  14  and  15  the  same  state  dis¬ 
tributions  are  shown  with  defining  spheres  and 
defining  ellipsoids.  Clearly  there  is  interference  in 
the  former  case,  but  no  interference  in  the  latter. 
Using  the  ellipsoids  to  define  the  state  volumes, 
we  performed  capacity  analyses  on  the  previous 
three  systems.  The  results  of  these  measures  are 


tabulated  in  Table  2,  where  we  can  see  that  the 
estimated  capacity  of  the  five-feature  silicon  sys¬ 
tem  has  nearly  doubled  to  having  27  of  32  states 
distinguishable. 


5.  Discussion 

Using  the  analyses  in  the  previous  sections,  we 
were  able  to  use  the  knowledge  from  the  three- 
feature  systems  to  establish  the  number  of  simula¬ 
tions  needed  for  higher  orders  of  features.  This  has 
allowed  us  to  reduce  the  simulation  times  neces¬ 
sary  for  investigating  other  systems.  Using  one  full 
set  of  simulations  for  a  five-feature  system,  we 
were  able  to  show  how  the  necessary  fraction  of 
neighbour  simulations  decreases  rapidly  as  the 
number  of  features  increases.  More  importantly, 
the  use  of  ellipsoids  rather  than  spheres  enables  a 
more  accurate  method  of  distinguishing  states, 
and  can  dramatically  improve  the  predictions  on 
capacity. 

To  compare  the  performance  of  this  alternative 
system  with  existing  technologies,  it  was  necessary 
to  estimate  the  achievable  system  density.  From 
the  simulations,  it  is  possible  to  calculate  a  linear 
density,  where  we  have  to  take  into  account  the 
spacing  of  250  nm  between  cells.  Assuming  that 
the  same  linear  density  can  be  achieved  in  the 
other  planar  dimension,  an  areal  density  can  be 
determined.  With  this  assumption,  the  maximum 
achieved  density  is  predicted  to  be  3.17  GBits/ 
cm2  for  the  silicon  five-feature  system.  Although 
this  is  much  less  than  our  system  goal,  it  is  still  five 
times  better  than  current  DVD  densities. 

Initial  simulations  with  other  degrees  of  free¬ 
dom  provide  some  evidence  indicating  even  higher 
densities  are  achievable.  By  allowing  each  feature 
to  take  on  more  than  two  possible  levels,  wc  can 
create  even  more  possible  states.  For  instance,  a 
five-feature,  four-level  system  would  have  1024 
possible  configurations  of  the  COI  (rather  than 
the  32  for  the  two-level  system).  A  larger  percent¬ 
age  of  states  will  be  indistinguishable,  but  it  could 
still  be  that,  for  instance,  128  are  distinguishable  - 
a  capacity  of  7  bits/A.  We  have  performed  a  very 
preliminary  analysis  of  such  a  system  -  using  the 
silicon  template  that  has  been  used  throughout 
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Fig.  13.  The  figure  on  the  left  shows  the  ploi  of  the  stale  vector  distributions,  for  every  COI*  for  the  silicon  three-feature  system,  along 
with  the  vector  indicating  the  centre  of  the  distribution.  The  figure  on  the  right  is  an  expanded  view,  without  the  centre  vectors, 
demonstrating  how  the  distributions  are.  in  general,  elongated  in  shape. 
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Fig.  14.  The  distribution  for  two  COTs  of  a  system,  and  the 
spheres  that  have  been  calculated  to  indicate  the  extent  of  the 
distribution.  Note  that  in  this  case,  to  determine  capacity,  both 
of  these  stales  would  be  discarded  because  of  the  overlap  in 
their  spheres. 


the  last  section  with  a  very  small  number  of 
neighbour  COIs.  Using  methods  that  over-estimate 
the  uncertainty  sphere  radii,  and  hence  provide  a 
conservative  estimate  of  the  capacity,  we  find  that 
at  least  54  of  the  1024  states  are  distinguishable. 
This  would  provide  roughly  double  the  linear 
density  of  the  best  system  discussed  so  far  in  this 
paper. 

Although  larger  densities  will  be  required  to 
surpass  current  magnetic  storage  technology,  there 
are  still  a  many  possibilities  to  explore  to  increase 
the  capacity  of  the  proposed  nanostructure  optical 
data  storage  method.  For  instance,  other  parame¬ 
ters  can  be  exploited  in  order  to  distinguish  states, 
such  as  the  wavelength  and  polarization  of  the 
source.  These  other  degrees  of  freedom  could  fea¬ 
sibly  dramatically  increase  the  capacity  of  the  sys¬ 
tem,  but  are  left  for  future  research. 

Besides  utilizing  methods  to  increase  the  de¬ 
grees  of  freedom,  there  are  many  avenues  to  be 
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Fig.  15.  The  same  two  distributions  as  in  Fig.  14,  except  with 
calculated  ellipsoids  rather  than  spheres  to  represent  their 
distributions.  In  the  capacity  measure  for  this  system,  these  two 
stales  would  now  be  valid  (distinguishable)  because  there  is  no 
overlap. 


Table  2 


A  comparison  of  the  capacity  of  the  three  different  systems, 
using  spheres  or  ellipsoids  to  define  the  state  distributions 


Capacity  (#  states//) 

Three- feature 

Five- feature 

Polymer 

Silicon 

Silicon 

Using  spheres 

3  of  8 

4  of  8 

13  of  32 

Using  ellipsoids 

6  of  8 

8  of  8 

27  of  32 

explored  for  improving  the  implementation  of 
the  system.  For  instance,  the  spatially  distributed 
near-field  power  data  could  be  converted  to 
spectral  data  for  measurement  in  the  far-field. 
One  could  imagine  utilizing  a  set  of  patterned 
array  of  organic  fluorophores  as  the  detector  le¬ 
vel  for  just  such  a  thing.  The  relative  spectral 
intensities  of  the  fluorophores  measured  in  the 
far-field  would  then  provide  the  data  measure 
for  a  COI. 


6.  Conclusion 

A  read-only  memory  system  was  proposed 
that  uses  near-field  optical  power  measurements 


to  distinguish  between  designed  deformities  in 
nanoscale  structured  cells.  Differing  deformities 
or  features  affect  the  power  measurements 
allowing  for  different  states  to  be  stored  and 
read.  The  capacity  (and  hence  the  density)  of  a 
given  system  is  governed  by  the  crosstalk  be¬ 
tween  cells,  causing  some  states  to  be 
indistinguishable. 

Measurements  and  methods  of  analysis  were 
developed  using  simple  simulations  of  three- 
feature  systems,  from  which  inferences  were  made 
on  the  fraction  of  simulations  necessary  at  higher 
orders  of  features.  Trends  showed  that  as  the  num¬ 
ber  of  features  increases,  the  fraction  of  simula¬ 
tions  necessary  to  get  an  accurate  estimate  of 
capacity  decreases.  In  the  example  of  a  five-feature 
system  that  was  presented,  approximately  20%  of 
the  full  set  of  simulations  was  enough  to  provide 
an  accurate  estimate  of  the  capacity.  Simulations 
for  three-feature  systems  also  provided  insight  on 
the  shape  of  the  power  distributions,  which  al¬ 
lowed  better  estimates  of  capacity  and  hence  den¬ 
sity.  Densities  of  3.17  GBits/cm2  have  been 
estimated  that  greatly  surpass  current  optical  stor¬ 
age  techniques. 

Although  the  densities  of  magnetic  storage  have 
not  yet  been  achieved,  the  optical  system  shows 
great  potential.  The  third  dimension  is  the  main 
factor  that  provides  hopes  for  larger  capacity  with 
an  increased  number  of  degrees  of  freedom,  such 
as  different  modes  of  polarization.  Another  possi¬ 
bility  is  investigating  the  effects  of  the  wavelength 
of  the  source  to  provide  additional  data  to  distin¬ 
guish  the  different  states.  Future  work  will  simu¬ 
late  the  system  in  three  dimensions  to  investigate 
these  additional  avenues  available  to  increase 
capacity. 
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Nonlinear  Codes  for  Dispersive  Nonlinear  Fibers 
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Abstract — We  consider  a  single  wavelength  channel  in  a  non¬ 
linear  dispersive  single-mode  fiber.  We  simulate  the  propagation 
of  binary  sequences  using  ON-OFF  keying  return- to-zero  modula¬ 
tion  at  a  bit  rate  of  10  Gb/s,  under  different  dispersion  and  input 
power  conditions.  We  determine  the  impact  that  dispersion  and 
nonlinearity  have  on  the  pulse  sequences  by  measuring  the  Eu¬ 
clidian  distances  (EDs)  among  all  pairs  of  received  sequences.  We 
then  model  the  set  of  EDs  as  a  fully  connected  graph  and  by  means 
of  a  clique-finding  algorithm  we  search  for  subsets  of  sequences 
that  form  codes  matched  to  this  nonlinear  channel,  resulting  in  sig¬ 
nificant  bit-error-rate  improvement. 

Index  Terms — Dispersive  nonlinear  fiber,  fiber-optic  communi¬ 
cations,  forward  error  correction,  nonlinear  codes. 


I.  Introduction 

G  ROUP-VELOCITY  dispersion  (GVD)  is  known  to  limit 
the  performance  of  fiber-optic  communication  systems 
by  causing  intersymbol  interference  (ISI)  that  increases  with 
propagation  distance  [I].  Moreover,  in  single-mode  fibers,  the 
pulse  propagation  is  power-dependent  due  to  Kerr  nonlinearity. 
This  nonlinearity  interacts  with  dispersion,  further  increasing 
the  pulsewidth  under  normal  dispersion  and  also  producing 
spectral  broadening  [2].  In  dispersion-compensated  systems, 
intrachannel  pulse  interaction  causes  distortion  and  jitter  that 
increases  with  power  and  propagation  length  [3]. 

In  this  letter,  we  present  a  numerical  analysis  of  the  effects 
that  GVD  and  nonlinearity  have  on  the  Euclidian  distances 
(EDs)  among  pulse  sequences  at  the  receiver  in  a  single  channel 
of  a  wavelcngth-division-multiplexing  (WDM)  system.  We  ob¬ 
serve  that  nonlinearity  and  dispersion  induce  a  significant 
decrease  in  the  minimum  ED  while  the  overall  ED  distributions 
are  not  considerably  changed.  We  describe  a  graph  theory-based 
methodology  with  which  we  obtain  codes  that  are  well  matched 
to  this  nonlinear  fiber  channel,  giving  a  significant  reduction 
in  bit-error  rate  (BER)  compared  to  linear  codes  of  the  same 
length. 


II.  DESCRIPTION  OF  THI-:  SIMULATIONS 

The  propagation  of  optical  pulses  is  governed  by  the  well- 
known  nonlinear  Schrodinger  (NLS)  equation  [4] 
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where  A  is  the  slowly  varying  amplitude  of  die  pulse,  z  is  the 
propagation  distance,  and  r  is  the  time  in  a  reference  frame 
that  travels  at  the  pulse  group  velocity.  GVD,  third-order  dis¬ 
persion,  nonlinearity,  and  losses  are  represented  by  the  parame¬ 
ters  /?2*/?3i7»  and  a>  respectively.  In  (l),  we  have  not  included 
higher  order  nonlinear  terms  like  those  that  account  for  stim¬ 
ulated  Raman  and  Brillouin  scattering,  as  they  are  negligible 
within  the  range  of  power  and  pulsewidth  we  consider  in  this 
letter.  To  solve  the  NLS  equation  given  in  ( 1),  we  use  the  sym¬ 
metrized  split-step  Fourier  method  [51,  [6]  with  constant  step 
size. 

We  simulate  the  propagation  of  10-bit  sequences  using 
Gaussian  pulses  and  ON-OE1  keying  retum-to-zero  modula¬ 
tion  at  a  bit  rate  of  10  Gb/s,  along  a  fiber  of  length  50  km. 
The  sequences  comprise  equally  .spaced  bit  slots,  with  each 
pulse  centered  in  its  bit  slot.  Also,  sequences  are  padded  with 
2.5  bit  slots  of  zeros  at  the  beginning  and  at  the  end  to  avoid 
energy  spilling  at  the  boundaries.  Sequences  are  necessary  to 
account  for  the  memory  in  the  channel,  manifested  through  the 
broadening  and  distortion  of  the  pulses,  which  can  extend  to 
several  adjacent  bit  slots.  We  find  that  using  10  bits  is  enough 
to  accurately  account  for  the  spreading  of  the  pulses,  while 
maintaining  a  reasonable  processing  time,  as  the  complexity 
of  our  procedure  grows  exponentially  with  the  number  of  bits. 
The  power  applied  to  the  pulses  ranges  from  1  to  100  mW. 

The  GVD  parameter  fa  is  varied  from  2  to  30  ps2/km,  in  the 
normal  dispersion  regime  and  we  fix  the  third-order  dispersion 
at  f%  =  0.1  ps:1/km.  The  nonlinear  parameter  7  is  set  to  a  fixed 
value  of  1.3  W  lknV_1.  We  assume  no  fiber  losses  in  order  to 
clearly  show  the  effects  of  nonlinearity  in  the  ED  measurements, 
although  we  do  consider  losses  when  evaluating  the  codes  as  we 
want  to  use  them  in  a  realistic  environment. 

In  order  to  analyze  the  impact  of  the  pulsewidth  for  a  given 
data  rate,  four  different  duty  cycles  are  studied,  namely.  2 57c. 
33%.  40%,  and  50%.  We  define  duty  cycle  as  the  full-width  at 
half-maximum  (FWHM)  of  the  pulse  intensity  over  the  bit  slot 
width.  For  a  given  peak  power  Po,  the  average  power  applied  to 
a  pulse  clearly  depends  on  the  duty  cycle.  Therefore,  to  make 
a  fair  comparison  of  performance  among  different  duly  cycles, 
we  apply  a  power  penalty  by  scaling  the  measured  ED  by  the 
average  intensity  of  each  pulse.  This  scaling  produces  a  nor¬ 
malized  ED  (NED)  that  also  allows  us  to  compare  the  effect  of 
nonlinearity  at  different  input  power  levels,  namely  1,  10.  and 
100  mW,  For  our  system,  we  assume  incoherent  and  direct  de¬ 
tection.  Hence,  ED  is  computed  in  the  intensity  domain 


where  IPt%  and  Iqj  are  the  measured  intensities  for  sequences  p 
and  q ,  respectively,  at  the  discrete  time  step  i.  The  summation 
is  over  all  samples  in  the  sequence. 
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Fig.  1.  (a)  Minimum  NED  versus  dispersion  ai  P(l  —  1, 10,  and  100  mW  for 

50%  duly  cycle.  (h>  Minimum  NED  versus  ai  Pa  —  100  mW  lor  several 
nonlinear  codes  using  a  50%  duly  cycle.  Hamming  (7,  4)  is  also  included  for 
comparison. 

Fig.  1(a)  shows  the  minimum  NED  among  pairs  of  sequences, 
versus  fa,  at  an  input  peak  power  Pq  of  1,  10,  and  100  mW, 
and  a  duty  cycle  of  50%.  Within  all  the  duty  cycles  considered 
in  our  simulations,  the  performance  at  50%  duty  cycle  consis- 
tently  showed  to  have  the  largest  minimum  NED.  We  observe 
in  Fig.  1(a)  that  the  minimum  NE!D  decreases  monotonically  as 
dispersion  increases.  As  the  power  is  increased,  the  minimum 
NED  decreases  more  rapidly  with  dispersion,  as  an  effect  of 
nonlinearity.  At  100  mW  and  30  ps2/km,  the  minimum  NED 
has  decreased  to  almost  a  half  of  the  normalized  value  found  at 
1  mW. 

111.  Nonlinear  Codes 

For  a  given  set  of  fiber  parameters  we  obtain  an  ED  map,  the 
set  of  EDs  among  all  pairs  of  10-bit  sequences.  Within  the  ED 
map,  we  search  for  subsets  of  sequences  with  large  minimum 
distance.  This  map  can  be  thought  of  as  an  interconnected  graph 
where  each  binary  sequence  is  a  node  and  every  ED  is  an  edge. 
This  graph  is  fully  connected,  because  there  is  a  nonzero  dis¬ 
tance  between  every  pair  of  nodes.  We  would  like  to  remove 
those  edges  with  an  ED  smaller  than  some  minimum  preferred 
value.  This  value  will  be  chosen  according  to  the  code  rate  we 
want.  Since  our  modulation  scheme  is  binary,  we  are  interested 
in  finding  codes  with  a  number  of  codewords  that  is  a  power 
of  two.  Since  it  is  not  possible  to  predict  the  minimum  ED  re¬ 
quired  to  produce  a  code  with  a  previously  chosen  number  of 
codewords,  we  have  to  do  an  iterative  search.  After  removing 
all  edges  with  ED  smaller  than  the  preferred  value,  we  seek  the 
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largest  subset  of  nodes  that  are  fully  connected,  i.e.,  where  there 
exits  an  edge  between  every  pair  of  nodes  in  the  subset.  This 
subset  of  nodes  forms  a  code.  This  is  known  as  the  clique  finding 
problem  [7],  [8],  and  only  an  exhaustive  search  can  guarantee 
the  optimal  result.  There  is  certain  symmetry  in  the  ED  map 
that  facilitates  the  search  for  a  good  (nonoptimal)  solution.  This 
solution  will  have  the  best  performance  under  the  fiber  parame¬ 
ters  in  which  it  was  found,  but,  as  we  will  show,  it  can  perform 
well  under  other  conditions. 

A  property  of  linear  codes  is  that  the  addition  of  any  two 
codewords  of  the  code  produces  a  codeword  that  is  part  of  the 
code.  A  nonlinear  code  is  a  more  general  set  of  codewords  that 
does  not  have  this  property.  As  a  consequence,  nonlinear  codes 
cannot  be  generated  by  a  matrix,  but  instead  by  a  lookup  table.  In 
practice,  codes  with  the  best  distance  properties  will  likely  to  be 
of  the  nonlinear  kind  due  to  the  nature  of  this  method,  and  also 
because,  under  1S1,  ED  does  not  necessarily  grow  monotonically 
with  binary  (Hamming)  distance. 

Fig.  1(b)  shows  the  minimum  NED  versus  fa  using  a  few 
of  these  nonlinear  codes,  for  Pq  —  100  mW  and  50%  duty 
cycle.  These  codes  are  10  bit  long  and  have  various  numbers 
of  codewords  according  to  their  code  rates  (the  rate  of  a  code 
is  R  =  log2  (#eode words/ 10).  They  are  all  optimized  for  fa  = 
30  ps2/km.  We  observe  that  there  is  a  significant  increase  of  the 
minimum  NED  for  the  entire  range  of  dispersion  compared  to 
the  minimum  NED  of  the  uncoded  sequences,  also  shown  in  the 
plot.  It  is  also  worthwhile  to  include  the  minimum  NED  per¬ 
formance  of  a  Hamming  (7.  4)  code,  shown  by  the  curve  with 
an  x  marker  in  Fig.  1(b),  as  this  shows  that  the  minimum  NED 
of  the  set  of  sequences  defined  by  this  code  decreases  signifi¬ 
cantly  as  dispersion  increases.  This  implies  that  the  code  does 
not  have  good  performance  under  LSI.  Table  1  shows  the  binary 
codewords  of  a  rate  0.6  nonlinear  axle  [which  we  label  NL1(10, 
6)),  optimized  for  P0  —  100  mW,  fa  =  30  ps2/km,  L  50  km, 
and  16.7  Gb/s  (such  that  the  information  rate  is  10  Gb/s).  Under 
these  fiber  conditions,  the  overall  minimum  NED,  the  maximum 
NED.  and  the  chosen  minimum  NED  to  get  the  code  NL1(1Q, 
6)  are  0.0526.  I ,  and  0.3409,  respectively  (all  values  normalized 
to  the  maximum  NED).  We  measure  the  BER  performance  of 
code  NLI(10,6)  in  Section  IV. 

iv.  BER  Perform ance 

Fig.  2  shows  the  BER  performance  of  the  axle  NLI  ( 10,  6)  at 
a  hit  rate  of  1 0  Gb/s,  with  parameters  fa  —  30  ps2/km,  fa  —  0. 1 
ps3/km,  Pq  =  100  mW.  a  =  0.2  dB/km,  7  —  1.3  W_1km_1, 
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Fig.  2.  BBR  performance  of  code  NU(I0,  6).  I  =  50  km,  P0  =  100  mW, 
G  —  30  ps'/km,  using  Gaussian  pulses  (hollow  markers)  and  hyperbolic  secant 
pulses  (tilled  markers). 
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Fig.  3.  BF’R  performance  of  code  Nl  1(  10,  6)  on  different  channel  conditions, 
using  Gaussian  pulses.  Po  =  20  mW,  L  —  100  km,  G  =  15,30.  and 
—30  ps-/km. 

and  L  =  50  km,  using  Gaussian  (hollow  markers)  and  hyper¬ 
bolic  secant  pulses  (filled  markers).  We  assume  additive  white 
Gaussian  noise  and  a  maxi  mum -likelihood  (ML)  receiver  with 
a  bandwidth  of  20  GHz.  The  performances  of  a  Hamming  (7, 
4)  code  and  that  of  the  uncoded  set  are  also  included.  A  power 
penalty  has  been  applied  to  each  curve  to  account  for  the  dif¬ 
ferent  code  rates,  and  the  coded  sequences  were  simulated  at 
16.7  Gb/s  so  that  the  information  bit  rate  was  maintained  at 
10  Gb/s.  We  see  in  Fig.  2  that  at  a  BER  =  10-7,  the  coding 
gain  of  this  code  is  about  5. 1  dB  over  the  uncoded  set  and  about 
3.2  dB  over  the  Hamming  (7,  4)  code,  which  has  a  comparable 
rate  of  0.57. 

Fig.  3  shows  the  BER  performance  of  the  code  NL1(10,  6) 
under  various  other  settings.  For  Pq  =  20  mW,  02  =  15  ps2/km, 
and  L  =  100  km  (bit  rate  and  other  parameters  kept  the  same 
as  in  previous  case),  the  coding  gain  is  1.5  dB  over  the  uneoded 
set  and  about  3.8  dB  over  the  Hamming  code,  at  a  BER  =  10"  ' 
(dashed  lines  with  hollow  markers).  Note  that  die  Hamming 
code  performance  is  worse  than  that  of  the  uncoded  sequences. 
This  occurs  because  the  Hamming  code  cannot  overcome  the  in¬ 
creased  LSI  that  the  pulses  suffer  as  they  propagate  at  a  higher  bit 


rate  (16.7  GB/s)  in  order  to  maintain  the  information  bit  rate  at 
10  Gb/s.  Under  the  same  fiber  parameters  but  with  a  dispersion 
of  30  ps2/km  (continuous  lines  with  black  markers),  the  results 
ofthe  code  NL  1(1 0,6)  are  excellent  too:  At  a  BER  =  10“  \  the 
coding  gain  is  5.8  dB  over  the  uncoded  bits,  and  3.8  dB  over  the 
Hamming  code. 

This  code  shows  to  give  good  results  in  the  anomalous  disper¬ 
sion  regime  as  well.  For  P0  =  20  mW  and  02  =  -30  ps2/km 
(maintaining  all  other  parameters  as  in  the  previous  ease),  the 
coding  gains  are  3.4  and  3.6  dB  over  the  Hamming  code  and 
the  uncoded  sequences,  respectively  (dotted  lines  with  gray 
markers). 

V.  Discussion  and  Conclusion 

Using  an  approach  from  graph  theory,  we  have  found  non¬ 
linear  codes  from  the  ED  maps  that  outperform  linear  codes 
of  similar  length  in  both  the  normal  and  anomalous  dispersion 
regime,  as  these  nonlinear  codes  account  for  the  nonlinearity 
and  dispersion  present  in  die  optical  fiber. 

We  have  determined  as  a  rule  of  thumb  that  good  performance 
can  be  achieved  with  code  NLI(10,  6)  if  die  ratio  of  disper¬ 
sion  length  Ln  =  ^fwhm/I/^I  over  propagation  length  L  is 
Ld/L  <  2.0  approximately.  Also,  an  ML  receiver  is  required. 
Other  codes  can  be  found  to  provide  good  performance  under 
lower  dispersion  and/or  bit  rates. 

These  codes  can  be  used  as  an  alternative  to  dispersion  com¬ 
pensation  and  can  be  used  in  WDM  systems.  This  mediodology 
could  still  be  used  on  systems  that  already  incorporate  disper¬ 
sion  compensation,  and  this  is  under  study. 

By  operating  at  high  power,  like  100  mW,  spectral  efficiency 
is  reduced  by  the  effects  of  nonlinearity,  as  it  causes  the  spec¬ 
trum  to  broaden.  In  practice,  spectral  broadening  is  not  con¬ 
spicuous  below  40  mW  in  the  presence  of  typical  fiber  losses 
and  within  the  dispersion  range  in  which  this  code  is  design 
to  give  good  performance.  Nevertheless,  code  NL1(  10,  6)  per¬ 
forms  very  well  at  lower  power,  so  it  is  not  essential  to  incur  in 
a  spectral  efficiency  penalty  to  obtain  good  results. 

We  are  working  on  increasing  the  code  length  by  extracting 
the  salient  features  of  the  nonlinear  codes  and  forming  linear 
codes  that  contain  these  features. 

References 

1 1 1  R.  J.  Nuyts,  Y.  K.  Park,  and  P.  Gal  lion,  “Performance  improvement  of 
10  Gb/s  standard  tiber  transmission  systems  by  using  the  SPM  elTect  in 
the  dispersion  compensating  liber.”  IEEE  Photon.  Technol.  Lett.,  vol.  K. 
pp.  1406-1408,  Oct.  1996. 

[2]  M.  Potasek,  G.  Agrawul,  and  S.  Pinault.  "Analytic  and  numerical  study 
of  pulse  broadening  in  nonlinear  dispersive  libers"!  Opt.  Sin  Amer.  H. 
vol.  3,  pp.  205-211.  1986. 

1 3 1  A.  Meeozzi,  C.  Clausen,  and  M.  Shtaif,  "Analysis  of  intrachannel 
nonlinear  effects  in  highly  dispersed  optical  pulse  transmission.”  IEEE 
Photon.  Technol.  Lett.,  vol,  12,  pp.  392-394,  Apr.  2000. 

1 4 1  G.  P.  Agrawal,  Nonlinear  Fiber  Optics.  San  Diego.  CA.  Academic, 
2001. 

151  - ,  Fiber-Optic  Communication  Systems.  New  York:  Wiley.  1997. 

[6J  T.  R.  Tuhu  and  M.  J.  Ablowilz,  “Analytical  and  numerical  aspects  of 
certain  nonlinear  evolution  equations,”  J.  Comput.  Pin  s .  vol.  55,  pp. 
203-230,  I9K4. 

1 7 1  M  R  Garey  and  D.  S.  Johnson.  “Complexity  of  near-optimal  graph  col 
oring.”  J.  ACM.  vol.  23,  no.  1.  pp.  43-19,  1976. 

[81  C.  lirun  and  J.  Kerbosh.  “Algorithm  457.  Finding  all  cliques  of  an  undi¬ 
rected  graph.”  Common.  ACM .  vol.  16,  no.  9,  pp.  575-577,  1973. 


Appendix  3 

Information  Propagation  in  a  Fast-Light  Medium 


NATURE  2016— 1 1/9/2003— VBICKNELL— 82343 


letters  to  nature 


Tlie  speed  of  information  in  a 
last-light’  optical  medium 

Michael  D.  Stenner  ,  Daniel  J.  Gauthier1  &  Marie  A.  Nerfeld 

1  Duke  University,  Department  of  Physics,  ami  The  Fitzpatrick  Center  for 
Photonics  and  Communication  Systems,  Durham,  North  Carolina  27708,  USA 
'Department  of  Electrical  and  Computer  Engineering,  The  Optical  Sciences 
Center,  University  of  Arizona,  Tucson,  Arizona  8572  L  USA 

One  consequence  of  the  special  theory  of  relativity  is  that  no 
signal  can  cause  an  effect  outside  the  source  light  cone,  the  space- 
time  surface  on  which  light  rays  emanate  from  the  source1. 
Violation  of  this  principle  of  relativistic  causality  leads  to  para¬ 
doxes,  such  as  that  of  an  effect  preceding  its  caused  Recent 
experiments  on  optical  pulse  propagation  in  so-called  ‘fast- 
light’  media — which  are  characterized  by  a  wave  group  velocity 
ug  exceeding  the  vacuum  speed  of  light  c  or  taking  on  negative 
values4 — have  led  to  renewed  debate  about  the  definition  of  the 
information  velocity  t  j.  One  view  is  that  i/j  =  i/g  (ref.  4),  which 
would  violate  causality,  while  another  is  that  vx  =  c  in  all 
situations9,  which  would  preserve  causality.  Here  we  find  that 
the  time  to  detect  information  propagating  through  a  fast-light 
medium  is  slightly  longer  than  the  time  required  to  detect 
the  same  information  travelling  through  a  vacuum,  even 
though  in  the  medium  vastly  exceeds  c.  Our  observations 
arc  therefore  consistent  with  relativistic  causality  and  help  to 
resolve  the  controversies  surrounding  superluminal  pulse 
propagation. 

The  speed  of  a  light  pulse  travelling  through  an  optical  material  is 
not  precisely  defined,  because  any  pulse  comprises  a  collection  of 
elementary  sinusoidal  waveforms,  each  with  a  distinct  frequency  o>. 
Each  constituent  sinusoid  travels  at  a  well-defined  velocity  known  as 
the  phase  velocity  up  —  c/n(a>),  where  n(o>)  is  the  refractive  index  of 
the  optical  material.  Approximate  theories  of  optical  pulse  propa¬ 
gation  predict  that  the  peak  travels  at  the  group  velocity 
Ug  —  c/(n  -f-  cad/i/dcaL^,,)  =  c/ng,  where  is  the  group  index 
and  o)0  is  the  central  frequency  of  the  wavepacket6. 

We  refer  to  the  quantity  dn/dca  as  the  dispersion  of  an  optical 
material.  For  typical  optical  materials,  there  exist  narrow  spectral 
regions  where  n(t*>)  is  a  decreasing  function  of  frequency  (that  is, 
dn/dco  <  0),  resulting  in  a  condition  known  as  anomalous  dis¬ 
persion  ;.  When  w0  is  within  such  a  region,  ng  can  be  less  than  one 
and  can  even  become  negative  when  the  anomalous  dispersion  is 
large.  This  results  in  ‘fast  light’,  for  which  it  is  possible  that  the  peak 
of  a  light  pulse  may  exit  the  optical  material  before  it  passes  through 
the  entrance  face14.  The  amount  of  fast-light  pulse  advancement  is 
largest  when  is  negative  and  near  zero  (Hg  large  and  negative). 

The  possibility  of  superluminal  group  velocities  (ug  >c  or 

<  0)  was  such  a  concern  to  researchers  around  1910  that  several 
conference  sessions  were  devoted  to  the  topic9.  Based  on  these 
discussions,  Sommerfeld  demonstrated  theoretically  that  the 
velocity  of  the  front  of  a  square  shaped  pulse  propagating  through 
any  medium  is  identically  equal  to  c  and  hence  relativistic  causality 
is  preserved10.  In  a  follow-up  study,  Brillouin  suggested  that  the 
group  velocity  is  not  physically  meaningful  when  the  dispersion  is 
anomalous  because  the  pulse  becomes  severely  distorted15.  More 
recent  research  investigating  the  propagation  of  smooth-shaped 
pulses  has  shown  that  this  conclusion  is  not  justified,  leading  to 
renewed  controversy*12-20. 

Another  outcome  of  the  discussions  in  the  early  1900s,  as 
recounted  in  the  preface  and  first  chapter  of  the  book  by  Brillouin9, 
was  a  reformulation  of  the  fundamental  postulate  of  the  special 
theory  of  relativity.  This  reformulation  states  that,  rather  than  limit¬ 
ing  the  speed  of  an  ‘object’,  it  is  the  information  velocity  v-x  that  is 


limited  by  c.  Unfortunately,  there  is  no  agreed-upon  definition  of 
the  information  velocity2. 

In  our  experiment,  we  use  a  fast-light  medium  that  exploits  the 
spectral  region  of  anomalous  dispersion  between  two  closely  spaced 
amplifying  resonances1^20  realized  by  creating  large  atomic  coher¬ 
ence21  in  a  laser-driven  potassium  vapour  as  shown  in  Fig.  la.  We 
obtain  larger  pulse  advancement  for  a  smooth  gaussian -shaped 
pulse,  as  shown  in  Fig.  lb,  in  comparison  to  the  experiment  of 
ref.  15,  by  increasing  the  gain  and  hence  the  size  of  the  anomalous 
dispersion.  The  larger  advancement  relative  to  the  pulse  width 
obtained  in  our  experiment  makes  it  easier  to  distinguish  the 
different  velocities  describing  pulse  propagation.  From  this  data, 
we  infer  that  =  — 19.6  ±  0.8,  indicating  that  we  are  operating  in 
the  highly  superluminal  regime. 

Measuring  u,  requires  an  understanding  of  the  fundamental 
mechanism  for  information  encoding  and  detection.  Garrison 
et  al 2  propose  that  new  information  is  encoded  on  an  optical 
pulse  by  creating  a  point  that  is  non-analytic  (for  example,  a 
discontinuity  in  the  pulse  amplitude  or  its  derivatives)  and  that 


Figure  1  Fast-light  pulse  propagation,  a.  Experimental  set-up.  The  potassium  vapours 
are  contained  in  two  uncoated  Pyrex  cells  ot  length  U2  =  20  cm  (to  suppress  unwanted 
parametric  instabilities^  and  heated  to  obtain  ari  atomic  number  density  of 
4.5  x  10M  atoms  cm  3  Linearly  polarized  coherence  preparation  laser  beams 
(frequencies  and  ut1 , )  are  combined  with  the  linear  and  orthogonally  polarized  pulses 
using  polarizing  beam  splitters.  The  pulses  are  detected  by  an  avalanche  photodiode  with 
a  25kHz  125  MHz  bandwidth.  The  coherence  preparation  beams  are  adjusted  with  u3d 
set  at  1 .36  GHz  to  the  high  frequency  side  of  the  centre  of  the  potassium  4 Pl/a  —  4Pf/a 
transition  and  -  u>a  -  23  MHz,  chosen  to  optimize  the  pulse  advancement  using 
procedures  similar  to  those  discussed  in  refs  1 5, 1 8  and  22.  The  pulses  are  generated  by 
passing  a  continuous  wave  laser  beam  through  an  acousto  optic  modulator  (AOM)  driven 
by  a  computer-controlled  arbitrary  waveform  generator.  The  time  origin  has  been  set 
arbitrarily  to  coincide  with  toe  peak  of  this  pulse,  b,  The  solid  line  shows  toe  temporal 
evolution  of  a  263,4-ns  long  (full-width  at  half-  maximum)  pulse  propagating  through  toe 
celts  when  the  lasers  are  tuned  far  from  the  atomic  resonance  and  hence  the  vapour-cell 
portion  of  the  path  is  equivalent  to  vacuum.  The  dashed  line  shows  the  observed  fast  lid  it 
pulse  advancement  for  a  smooth  pulse  shape  when  the  coherence- preparation  laser 
is  tuned  near  the  atomic  resonance  and  o>  0  is  set  between  the  gain  resonances.  The 
peak  of  toe  pulse  is  advanced  by  Gdv  =  27.4  ns  ±  1.1  ns,  corresponding  to  a  relative 
pulse  advancement  of  10.4%.  Using  t **  -  L/c  Llv^  with  Uc  -  1.3  ns,  we  find 
V'jc  =  -0.051  t  0.002.  Careful  inspection  of  toe  fast-light  pulse  reveals  that  it  has 
been  compressed  by  1 .9%,  which  is  due  primarily  to  the  frequency  dependence  of  the 
gain 181 
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this  point  always  travels  at  c  regardless  of  the  value  of  the 
other  velocities  associated  with  the  pulse  3.  Essentially,  they  have 
generalized  Sommerfeld’s  concept  of  the  front  velocity  to  a  non- 
analytic  point  of  the  pulse  amplitude,  where  the  front  of  a  square¬ 
shaped  pulse  is  an  example  of  a  point  of  non  analyticity.  These 
workers*'  ■  suggest  that  the  point  of  non-analyticity  is  the  only  part 
of  the  pulse  representing  new  information  because  measurements  of 
the  early  part  of  the  pulse  cannot  be  used  to  predict  anything  about 
the  part  of  the  pulse  arriving  after  the  point  of  non-analyticity,  and 
hence  v ,  equals  the  speed  of  a  point  of  non-analyticity.  For  counter¬ 
views,  see  refs  2,  4,  23  and  24.  We  note  that  some  aspects  of  their 
proposal  have  been  verified  using  electronic  circuits  where  no 
propagating  waves  are  involved,  and  hence  only  issues  of  causality, 
but  not  relativistic  causality,  can  be  tested25*-27. 

To  enhance  our  ability  to  estimate  the  location  of  this  non- 
analytic  point  in  the  presence  of  noise,  we  use  two  optical  pulses  that 
are  initially  identically  gaussian-shaped,  which  allows  us  smoothly 
to  turn  on  the  pulse  amplitude  to  a  level  above  the  noise  floor  of  our 
detection  electronics  and  to  monitor  the  fast -light  pulse  advance¬ 
ment.  Near  the  peak  of  the  gaussian  function  and  at  the  same 
moment  for  both  symbols,  we  switch  the  amplitude  of  the  gaussian 
function  to  a  high  (i)  or  low  (0)  value  for  the  remainder  of  the 
pulse.  The  moment  when  a  decision  is  made  to  switch  between  the 
symbols  corresponds  to  the  point  of  non-analyticity.  Note  that  this 
transition  is  smoothed  out  by  the  finite  response  time  of  the  optical 
switch. 

The  location  of  the  point  of  non-analyticity  is  determined  by 
detecting  the  arrival  of  new  information  using  a  receiver  that  can 
distinguish  between  symbols  to  a  desired  level  of  certainty,  charac¬ 
terized  by  the  bit  error  rate  (BER).  Before  the  arrival  of  the  point  of 
non-analyticity  at  the  detector,  we  expect  no  detected  information, 
corresponding  to  a  BER  of  1/2.  Once  the  point  of  non-analyticity 
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Figure  2  Transmitting  information -encoded  optical  pulses  through  a  fast-light  medium, 
a,  Transmitting  '0'  and  '1  *  through  the  fast-light  medium  (dashed  line)  and  vacuum  (solid 
line).  Each  symbol  is  transmitted  separately  through  the  medium  and  vacuum,  where  each 
curve  is  an  average  of  50  pulses,  b,  High -resolution  plot  of  part  of  a  with  an  offset  for 
clarity.  The  amplitude  of  the  advanced  and  vacuum  pulses  have  been  scaled  so  that  their 
heights  would  be  the  same  if  a  gaussian  pulse  propagated  through  the  system,  as  in 
Fig.  1  b.  The  error  bar  indicates  the  typical  standard  deviation  of  the  pulse  amplitudes. 
From  Fig.  2a,  it  is  seen  that  the  fast- light  medium  advances  the  early  part  of  the  pulses 
during  the  smooth  turn  on,  identically  to  that  observed  for  the  full  gaussian  -shaped  pulses 
shown  in  Fig.  1b.  Most  important  is  the  observation  that  both  symbols  are  the  same  for 
early  times  so  that  it  is  not  possible  to  distinguish  between  them.  Hence,  no  information 
can  yet  be  conveyed  to  a  receiving  parly  at  the  end  of  the  communication  channel. 


propagates  past  the  detector,  the  received  information  will  grow 
smoothly  from  zero  and  the  BER  drops.  A  symbol  is  considered  to 
be  detected  when  the  BER  falls  below  some  threshold.  Hence,  the 
detection  time  of  information  is  later  than  the  time  when  infor¬ 
mation  is  first  available  at  the  detector,  even  for  pulses  propagating 
through  vacuum.  This  detection  latency  At  depends  on  the  charac 
teristics  of  the  medium  through  which  the  pulses  propagate,  the 
shape  of  the  symbols,  the  detection  algorithm,  noise*’ 1  in  the 
detection  process  (including  quantum  noise'-  "’h  and  the  BER 
threshold.  It  increases  as  the  signal -to -noise  ratio  decreases  because 
it  takes  longer  for  the  receiver  to  achieve  the  same  BER.  Achieving 
the  limit  At  — ■ ►  0  requires  the  use  of  optimal  symbol  shapes  and 
detection  algorithms,  and  infinite  energy  in  the  optical  pulse  so  that 
the  signal-to-noise  ratio  of  the  detected  waveform  approaches 
infinity.  Although  it  is  possible  to  estimate  Af  for  a  specific 
experimental  apparatus,  it  cannot  be  measured  directly  because  it 
requires  measuring  the  point  of  non-analyticity.  A  crucial  aspect  of 
our  experiment  is  to  make  Af  as  small  as  possible  and  to  make  it  as 
similar  as  possible  for  both  vacuum  and  advanced  pulses. 

Figure  2  shows  the  propagation  of  both  symbols  through  the  fast- 
light  medium  and  vacuum.  From  a  simple  visual  inspection  of  the 


Figure  3  Detecting  the  arrival  of  new  information.  Shown  is  the  BER  as  a  function  of  tin 
upper  limit  of  the  integration  time  r  for  the  vacuum  (solid  line)  and  advanced  (dashed  line) 
pulses.  The  horizontal  line  indicates  the  detection  threshold,  a,  Experimental 
observations.  The  time  origin  has  been  selected  arbitrarily,  b,  Theoretical  predictions 
based  on  Maxwell's  equations.  The  time  origin  corresponds  to  the  moment  when  the  point 
of  non-analyticity  first  arrives  at  the  detector.  The  BER  is  determined  using  the  following 
matched- filter  procedure.  For  each  symbol  (0  and  1),  there  are  50  pulse  waveforms. 
Reference  waveforms  for  each  symbol,  denoted  by  Rn(t)  and  R  ,(f).  are  generated  by 
averaging  49  pulse  waveforms  for  each  symbol  and  the  integral  0(7)  -  /0(r)  -  /,  (r)  is 
determined  for  each  pulse,  where  the  waveform  being  detected  is  not  included  in  the 
respective  reference  waveform.  Here:  /,(r)  =  Jfr; f  r  x(  t)R,(  t  )6t/a,  (r(f  )A//(  t)  (for  /  0, 

1),  where  40  is  an  individual  pulse  waveform,  f.;  is  an  integration  start  lime  chosen 
arbitrarily  at  a  time  in  the  early  beginning  of  the  pulse,  N,(t)  -  J  r'  '  fljV  )df. 
a;(rrt)  =  Jfr,+Ta  and  t„  is  a  normalization  integration  time  chosen 

arbitrarily  on  the  leading  edge  of  the  gaussian  pulse  before  the  point  where  the  symbols 
separate.  For  each  symbol  type,  the  probability  density  of  D  is  fitted  to  a  gaussian 
dtetribution  normalized  to  an  area  of  1/2  and  the  overlapping  area  of  these  two  gaussian 
distributions  is  the  BER. 
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data,  we  see  that  the  time  where  it  is  possible  to  first  distinguish 
between  the  two  symbols  for  the  advanced  pulses  is  slightly  later 
than  the  time  where  it  is  possible  to  first  distinguish  between  the 
same  two  symbols  propagating  through  vacuum.  In  addition,  it  is 
seen  that  the  manner  in  which  the  average  symbol  waveforms 
separate  for  the  vacuum  and  advanced  case  are  only  slightly 
different,  so  that  the  detection  latency  times  (denoted  by  Afv.ac 
and  Afjdv  respectively)  should  be  similar. 

To  quantify  our  results,  we  determine  the  BER  for  the  vacuum 
(Fig.  3a,  solid  line)  and  advanced  (dashed  line)  pulse  pairs  using  an 
integrate-and-dump  matched  filter  technique.  The  BER  is  high  for 
final  observation  times  in  the  range  between  —40  and  —25  ns, 
during  which  the  pulse  amplitudes  are  large  (see  Fig.  2a).  Hence, 
even  though  the  $igna!-to-noi$e  ratio  for  a  single  pulse  is  high  at 
these  times,  the  pulses  are  not  yet  distinguishable  and  no  infor¬ 
mation  is  detected.  Placing  the  detection  threshold  at  BER  =  0.1, 
chosen  to  keep  Afvmc  and  Arddv  small,  we  determine  the  detection 
time  for  vacuum  (advanced)  pulse  pairs  Tvac  (7* dv)  and  the 
difference  in  detection  times  7',  —  T,jv  —  7*ViK.  The  time  difference 
is  approximately  constant  for  BER  values  around  0.1;  its  average 
value  in  the  range  of  BERs  between  0.08  and  0.2  is  equal  to 
3.2  ±  1.5  ns.  Our  observations  demonstrate  that  the  information 
detection  time  for  pulses  propagating  through  the  fast-light  med¬ 
ium  is  longer  than  the  detection  time  for  the  same  information 
propagating  through  vacuum,  even  though  the  group  velocity  is  in 
the  highly  superluminal  regime  for  the  fast-light  medium. 

From  our  direct  measurement  of  Tj,  we  do  not  know  whether  an 
observed  difference  between  the  detection  times  is  due  to  changes  in 
the  detection  latencies  or  differences  in  the  information  velocities 
for  vacuum  (u,iVttC)  and  the  fast-light  medium  (uuadv)-  The  relation 
among  these  quantities  is  given  by: 

7  l  =  (E/u^julv  f*/ ^i.vac)  T  (AfjJv  Arvac)  (0 

To  gain  some  insight  about  the  importance  of  detection  latency, 
we  analyse  a  mathematical  model,  based  on  Maxwell’s  equations, 
that  describes  approximately  the  generation,  propagation,  and 
detection  of  our  symbols.  Consistent  with  previous  research2*5, 
this  model  predicts  that  U|,aav  ~  Ui.vac  ~  c>  and  hence  7]  is  com 
pletely  determined  by  (A fadv  —  Afvaic).  Using  the  same  matched- 
filtering  approach,  we  determine  the  predicted  BER  as  shown  in 
Fig.  3b.  We  see  that  information  is  detected  later  for  the  advanced 
pulses  than  for  the  vacuum  pulses,  qualitatively  similar  to  the 
experimental  observations.  We  find  that  7’j  =  1.5  ±  0.5  ns,  where 
the  error  only  accounts  for  statistical  uncertainty  in  the  BER 
determination.  The  fact  that  Tty*  0  demonstrates  that  subtle 
changes  in  the  shape  of  the  symbols  after  information  has  been 
encoded  give  rise  to  substantial  changes  in  the  detection  latency.  The 
predicted  value  of  T,  might  be  smaller  than  the  observed  time  owing 
to  our  assumption  that  the  fast-light  medium  does  not  change  the 
noise  properties  of  the  optical  pulses2*. 

Using  the  model  prediction  for  (A ftldv  —  Af^-)  in  equation  (I) 
and  taking  v ,(V,C  —  e,  we  find  that  i/Uav  “  (0.4  ±  0.5)c.  Thus,  our 
observations  arc  consistent  with  the  special  theory  of  relativity  even 
for  a  medium  where  ug  is  highly  superluminal,  demonstrating  that 
the  peak  of  the  advanced  pulse  at  the  exit  face  of  the  medium  (see 
Fig.  Ib)  is  not  causally  connected  to  the  peak  at  the  entrance  face2. 


Because  our  analysis  makes  no  assumptions  about  the  sources  of 
noise  in  the  encoding,  transmission  and  decoding  process,  our 
general  experimental  approach  and  conclusions  should  hold  even  in 
the  limit  where  quantum  fluctuations  are  dominant □ 
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Information  Propagation  in  a  Slow-Light  Medium 


Fast  (but  causal)  information  trasmission  in  an  optical 
material  with  a  slow  group  velocity 
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We  demonstrate  transmission  of  information  encoded  on  optical 
pulses  at  a  velocity  60  times  faster  than  the  group  velocity  which 
describes  the  speed  of  the  peak  of  the  pulses.  Our  experiments 
are  conducted  using  a  laser-driven  potassium  vapor  with  vg  ~  0.01c, 
where  c  is  the  speed  of  light  in  vacuum.  We  measure  the  veloc¬ 
ity  of  information  v%  by  transmitting  different  pulse  shapes  through 
the  vapor  and  measuring  the  time  required  to  distinguish  the  dif¬ 
ferent  shapes.  Our  observation  that  vg  directly  contradicts  the 
conventional  wisdom  that  u*  =  vg  when  0  <  vg  <  c. 

One  of  the  most  fundamental  characteristics  of  a  pulse  of  electromagnetic'  radiation 
is  the  speed  at  which  it  propagates.  In  vacuum,  the  pulse  travels  at  speed  c  and  does 
not  change  its  shape.  On  the  other  hand,  a  pulse  traveling  through  a  dispersive  material 
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(such  as  an  optical  fiber  or  a  gas  of  atoms)  experiences  some  degree  of  pulse  reshaping 
due  to  the  frequency  dependence  of  the  refractive  index  n(u).  The  reshaping  arises  from 
the  fact  that  a  pulse  is  a  coherent  superposition  of  component  sinusoidal  waves  (/),  each 
of  which  travels  through  the  medium  at  a  different  speed  known  as  the  phase  velocity 
vv  —  c/n(w). 

In  general,  there  is  no  simple  way  to  completely  characterize  the  speed  of  a  pulse  duo 
to  the  reshaping  effect.  To  address  this  problem,  it  is  customary  to  introduce*,  several 
characteristic  velocities  that  describe  different  aspects  of  the  reshaping  process.  For  ex¬ 
ample,  it  is  well-known  that  the  peak  of  the  pulse  travels  approximately  at  the  group 
velocity  vg  =  c/(n  +  wdn/duj\u=u)0),  where  uq  is  the  central  frequency  of  the  pulse.  (1,2). 
This  velocity  describes  the  lowest-older  effect  of  the  reshaping  process:  a  translation  in 
space-time  of  the  pulse  in  comparison  to  the  same  pulse  traveling  in  vacuum. 

Recent,  research  has  shown  that  it.  is  possible  to  make  dramatic  changes  in  the  group 
velocity  of  a  gas  of  atoms  by  driving  it  with  an  intense  laser  beam,  thereby  creating 
large  atomic  coherence,  which  in  turn  induces  large  changes  in  cln/du  (.7).  For  example, 
by  inducing  large  normal  dispersion  (positive  dn/cLj)  or  anomalous  dispersion  (negative 
dri/du) ),  it  is  possible  to  slow  light  to  less  than  20  m/s  (4),  stop  or  freeze,  light  (5  7),  or 
even  make  superluininal  pulses  (6\  9). 

Superluminal  pulses,  with  vg  >  c  or  vg  <  0,  have  created  considerable  controversy 
because  their  existence  might  suggest  that  information  can  travel  faster  than  c,  which 
would  violate  relativistic  causality  (10).  This  controversy  has  been  addressed  in  a  recent 
experiment  where  it  was  found  that  the  information  velocity  in  fast-light  media  is  indeed 
less  than  c  (11). 

In  contrast ,  there  has  been  little  research  to  directly  measure  the  information  veloc¬ 
ity  in  a  normally-dispersive,  slow-light  medium.  This  situation  is  perhaps  understand- 
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able*:  many  textbooks  on  electromagnetism  demonstrate  theoretically  that,  in  a  normally- 
dispersive  material,  vg  is  approximately  equal  to  the  signal  velocity  v3  as  defined  by  Som- 
nu'jfeld  and  Drillouin  nearly  a  century  ago  (1.2).  The  name  “signal  velocity'’  implies  that 
it  is  related  to  the  transmission  of  information,  a  concept  that  is  echoed  in  many  texts 
and  research  articles  (12.  IS).  As  a  result.,  it  is  commonly  believed  that  vt  —  vs  as  vg  for 
normally-disporsive  media. 

Cliiao  and  collaborators  offer  a  different  view.  They  show  that  va  can  exceed  c, 
suggesting  that  va  cannot  be  the  true  velocity  of  information  (10).  Also,  the  intensity 
threshold  used  to  determine  va  is  arbitrary  (2),  allowing  for  different  values  of  vs  for  the. 
same  pulse.  As  an  alternative,  they  propose  that  information  is  contained  only  in  points 
of  non-analytieity  on  electromagnetic  waveforms,  and  therefore  that  vt  is  equal  t.o  the 
speed  with  which  these  points  propagate  (10, 14).  An  example  of  such  a  lion-analytic 
point  is  a  discontinuity  in  the  waveform  or  one  of  its  derivatives.  Their  proposal  is  based 
on  the.  notion  that  measurements  of  a  waveform  on  one  side  of  a  non-analytic  point  cannot 
be  used  to  infer  anything  about  the  waveform  on  the  other  side  of  the  point.  Hence,  the 
arrival  of  a  noil-analytic  point  brings  new  information.  Perhaps  the  most,  elegant  aspect 
of  their  proposal  is  that  it  can  be  viewed  as  an  extension  of  Sommerfold  and  Brillouin's 
earlier  research  where  they  showed  that  the  front  velocity  vj  (the  velocity  at  which  the 
lending  edge  of  a  waveform  propagates)  is  always  equal  to  c  regardless  of  the  properties  of 
the  medium  (:J).  In  linear  media,  one  can  always  view  a  waveform  with  multiple  points  of 
non-analvl  icit.y  as  a  superposition  of  multiple  waveforms  that  me  each  analytic  after  an 
initial  front.  Therefore,  according  to  the  theory  of  Cliiao  and  collaborators,  vt  =  Vf  =  c. 

In  order  to  help  resolve  the  conflict  between  Chiao’s  new  theory  and  the.  common  belief 
that  vs  -  we.  observe  the  propagation  of  information  in  a  medium  where  vs  &  vg  <  c. 
We  prepare  this  “slow-light”  medium  by  creating  a  single  amplifying  resonance,  realized 
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by  creating  large  atomic  coherence  in  a  laser  driven  potassium  vapor  (Fig.  la).  In  the 
spectral  region  of  the  amplifying  resonance,  the  group  velocity  is  very  slow,  and  pulses 
injected  with  center  frequencies  in  this  region  experience  large  delays  in  comparison  to 
identical  pulses  traveling  through  vacuum.  Figure  lb  shows  pulses  propagating  through 
vacuum  (solid  line)  and  through  the  slow-light,  medium  (dashed  line).  We  find  that  the 
peak  of  the  pulse  is  delayed  by  t^i  =  67.5  ±  2  ns,  corresponding  to  a  relative  pulse  delay 
of  25%.  Using  t^i  =  L/vg  —  L/c  with  L/c  =  0.66  ns.  we  infer  that  r0/c  —  0.01  ±  0.0003. 

To  enhance  our  ability  to  estimate  tin1  location  of  a  point  of  non-analyticity  in  the 
presence  of  noise,  we  use  two  optical  pulse  shapes  (our  “symbols')  that  have  identically 
Gaussian-shapocl  leading  edges.  This  allows  us  to  smoothly  turn  on  the  pulse  amplitude 
to  a  level  above  the  noise  floor  of  our  detection  electronics  and  to  monitor  the  pulse  delay. 
Near  the  peak  of  the  Gaussian  function  and  at  the  same  moment  for  both  symbols,  we 
switch  the  amplitude  of  the  Gaussian  function  to  a  high  ("1")  or  low  (“0”)  value  for  tin' 
remainder  of  the  pulse.  The  moment  when  a  decision  is  made  to  switch  between  the 
symbols  corresponds  to  the  point  of  non-analyticity.  Note  that  this  point  is  smoothed  out 
by  the  finite  rise,  time  of  the  optical  switch. 

Figure  2  shows  the  propagation  of  both  symbols  through  the  slow-light  medium  and 
vacuum  (2a)  along  with  an  enlargement  (2b)  of  the  diagram  in  tin1  vicinity  of  the  transition 
between  the  symbols.  From  Fig.  2a.  it  is  seen  that  the  slow  light  medium  delays  the  early 
part  of  the  pulses  during  the.  smooth  turn  on,  identically  to  that  observed  for  the  full 
Gaussian-shaped  pulses  shown  in  Fig.  lb.  Most  important  is  the  observation  that  both 
symbols  are  the  same  for  early  times  so  that  it  is  not  possible  to  distinguish  between 
them.  Hence,  no  information  can  yet,  be  conveyed  to  a  receiving  party  at  the  end  of  t  he 
communication  channel.  From  a  simple  visual  inspection  of  the  data,  we  see  that  the  time 
when  it  is  possible  to  first  distinguish  between  the  tw'o  symbols  for  the  delayed  pulses  is 
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nearly  the  same  as  the  time  when  it  is  possible  to  first  distinguish  between  the  same  two 
symbols  propagating  through  vacuum. 

The  arrival  time  of  the  information  is  determined  by  observing  the  pulses  with  a 
receiver  that  can  distinguish  between  the  two  symbols  to  a  particular  level  of  certainty, 
characterized  by  the  bit -error-rate  (BER).  Before  the  arrival  of  the  point  of  non-analvticity 
at  the  detector,  we  expect  no  detected  information,  corresponding  to  a  BER  of  1/2.  Once 
the  point  of  nou-analyticity  propagates  past  the  detector,  the  BER.  drops  as  the  received 
information  grows  smoothly  from  zero.  A  symbol  is  considered  to  be  detected  when  the 
BER  falls  below  some  threshold.  Hence,  the  detection  time  of  information  is  later  than 
the  time  when  information  is  first  available  at  the  detector,  even  for  pulses  propagating 
through  vacuum.  This  detection  latency  At  depends  on  the  characteristics  of  the  medium 
through  which  the  pulses  propagate,  the  shape  of  the  symbols,  the  detection  algorithm, 
noise  in  the  detection  process  (including  quantum  noise  (15.16)),  and  the  BER  threshold. 
It  increases  as  the  signal- to- noise  ratio  decreases  because  it  takes  longer  for  the  receiver 
to  achieve  the  same  BER  .  Achieving  the  limit  A t,  — >  0  requires  the  use  of  optimal  symbol 
shapes  and  detection  algorithms,  and  infinite  energy  in  the  optical  pulse  so  that  the 
signal-to-noise  rat  io  of  the.  detected  waveform  approaches  infinity.  Although  it  is  possible 
to  estimate  At  for  a  specific  experimental  apparatus,  it  cannot  be  measured  directly 
because  it  requires  measuring  the  point  of  non-analytirity.  Making  At  very  small  and 
very  similar  for  both  vacuum  and  delayed  pulses  is  a  primary  consideration  in  our  choice 
of  symbol  shapes.  In  Fig.  3b,  we  sec  that  the  manner  in  which  the  average  symbol 
waveforms  separate  for  the  vacuum  and  delayed  case  are  only  slightly  different  so  that 
the  detection  latency  times  (denoted  by  A tvac  and  A respectively)  should  be.  similar. 

To  quantity  our  results,  we  determine  the  BER  for  the  vacuum  (Fig.  2a.  solid  line)  and 
delayed  (dashed  line)  pulse  pairs  using  an  integrate-and-dump  matched  filter  technique 
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(17).  It  is  seen  that  the  BER  is  high  for  final  observation  times  in  the  range  between 
—40  and  —20  ns  during  which  the  pulse  amplitudes  are  large  but  not  yet  distinct  within 
the  noise  limit  (see  Fig.  2a).  Even  though  the  signal-tonoise  ratio  for  a  single  pulse  is 
high  at  these  times,  the  pulses  are  not  yet  distinguishable  and  therefore  no  information  is 
detected.  Placing  the  detection  threshold  at  BER  =  0.1,  chosen  to  keep  A/,,ac  and  At*/ 
small,  we  determine  the  detection  time  for  vacuum  (delayed)  pulse  pairs  Tvac  ( 7*/ )  and 
the  difference  in  detection  times  T,  —  Tdel  —  Tvac.  The  time  difference  is  approximately 
constant  for  BER  values  around  0.1;  its  average  value  in  the  range  of  BERs  between  0.08 
and  0.2  is  8.5  ±  0.5  ns. 

Based  only  on  a  direct  measurement  of  T*,  one  cannot  determine,  whether  an  observed 
difference  between  the  detection  times  is  due  to  changes  in  the  detection  latencies  or 
differences  in  the  information  velocities  for  vacuum  (vliVac)  and  the  slow-light,  medium 
(vi,dei)-  The  relation  among  these  quantities  is  given  by 

1\  =  (L/Videl  L/Vi%vac)  T  (A tdei  A tvac).  (!) 

To  gain  some  insight  about  the  importance  of  detection  latency,  we  analyze  a  mathemat¬ 
ical  model,  based  on  Maxwell’s  equations,  that  describes  approximately  the  generation, 
propagation,  and  detection  of  our  symbols.  Consistent  with  previous  research,  this  model 
predicts  that  ry*/  —  vliVac  =  c  (2,  10,  12),  and  hence  T,  is  completely  determined  by 
(At*/  -  A tvac).  Using  the  same  matched  filtering  approach,  we  determine  the  predicted 
BER  as  shown  in  Fig.  3b.  We  see  that  information  is  detected  later  for  the  delayed 
pulses  than  for  the  vacuum  pulses,  qualitatively  similar  to  the  experimental  observations. 
We  find  that  T»  -  8  ±  2  ns,  where  the  error  arises  from  uncertainty  in  the  similarity 
between  the  model  and  physical  experiment.  The  fact  that  T,  7^  0  demonstrates  that 
subtle  changes  in  the  shape  of  the  symbols  after  information  has  been  encoded  give  rise 
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to  substantial  changes  in  the  detection  latency. 

Using  the  model  prediction  for  (A t^i  —  Atmc)  in  Eq.  1  and  taking  vlA,ac  =  c,  our 
best  estimate  (18)  is  that  Vi&i  —  0.6c,  which  is  60  times  faster  than  vg.  Thus,  even 
in  normally-* lispersive  slow-light  media,  information  does  not  propagate  at  the  group 
velocity.  This  result  is  in  direct  contradiction  to  the  conventional  wisdom  regarding 
propagation  of  information  on  optical  pulses  and  hopefully  helps  resolve  a  century-old 
confusion.  Additionally,  this  result  serves  as  a  demonstration  of  a  new  technique  for 
measuring  Vi.  Future  refinements  of  this  technique  will  allow  precise  examination  of 
Chiao’s  proposal  that  the  information  velocity  is  always  equal  to  c. 
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Materials  and  Methods 
Fig.  SI 


Fig.  1.  Slow  light,  pulse  propagation,  (a)  Experimental  setup  (17).  (b)  Temporal 

evolution  of  265-ns-long  (full  width  at  half  maximum)  Gaussian-shaped  pulses  propagating 
through  the  slow  light  medium  (dashed  line)  and  vacuum  (solid  line),  demonstrating  slow- 
light  pulse  delay.  The  time  origin  has  been  set  arbitrarily  to  coincide  with  the  peak  of  the 
vacuum  pulse.  Careful  inspection  of  the  delayed  pulse*  reveals  that  it,  has  been  expanded 
to  a  pulse  width  of  311  ns  (17%  expansion),  which  is  due  primarily  to  the  frequency 
dependence  of  the  gain. 

Fig.  2.  (a)  Transmitting  “0”  and  “1”  through  the  slow  light  medium  (dashed  line) 
and  vacuum  (solid  line).  Each  symbol  is  transmitted  separately  through  the  medium 
and  vacuum,  and  each  curve  is  an  average  of  50  pulses,  (b)  High-resolution  plot  of  (a). 
The  amplitude  of  the  advanced  and  vacuum  pulses  have  been  scaled  so  that  their  heights 
would  be  the  same  if  a  Gaussian  pulse  propagated  through  the  system,  as  in  Fig.  lb.  The 
error  bar  indicates  the  typical  standard  deviation  of  the  pulse  amplitudes. 

Fig.  3.  Detecting  the  arrival  of  new  information:  the  BER  as  a  function  of  the  upper  limit, 
of  the  integration  time  for  the  vacuum  (solid  line)  and  advanced  (dashed  line)  pulses.  The 
horizontal  line  indicates  the  detection  threshold,  (a)  Experimental  observations.  The  time 
origin  is  the  same  as  in  Fig.  lb.  (b)  Theoretical  predictions  based  on  Maxwell’s  equations. 
The  time  origin  corresponds  to  the  moment  when  the  point,  of  non-ana  Ivtieity  first  arrives 
at  the  detector. 
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Supporting  Online  Material 


Preparation  of  the  slow-light  medium 

Tlu'  potassium  vapor  is  contained  in  an  uncoated  pyrex  cell  of  length  L  20  cm 
and  heated  to  obtain  an  atomic  number  density  of  3.5  x  1012  atoms/ cm3.  A  linearly- 
polarized  coherence- preparation  laser  beam  (frequency  u>j)  is  combined  with  the  linear 
and  orthogonally-polarized  pulses  (center  frequency  u/q  fa  uij  +  462  MHz)  using  a  polar¬ 
izing  beam  splitter.  The  pulses  are  generated  by  passing  a  continuous-wave  laser  beam 
through  an  acousto-optic  modulator  (AOM)  driven  by  a  computer-controlled  arbitrary 
waveform  generator.  The  pulses  are  detected  by  a  photoreceiver  with  a  25  kHz  -  125  MHz 
bandwidth. 

The  vacuum  pulses  are  observed  by  tuning  the  laser  very  far  from  the  atomic  resonance 
so  that,  neither  the  coherence-preparation  laser  beam  nor  the  pulses  interact  with  the 
atoms.  In  this  configuration,  the  potassium  vapor  is  equivalent,  to  vacuum. 

To  observe  slow  light  pulse  propagation,  the  coherence  preparation  beam  is  adjusted 
with  u>d  set  at  2.61  GHz  to  the  high-frequency  side  of  the  center  of  the  potassium  4S'I/2  — » 
4P[/2  transition.  This  creates  an  amplifying  resonance  at.  frequency  w,/  +  462  MHz,  the 
ground-state  splitting  of 39 ITS.  We  then  set  the  pulse  frequency  to  the  center  frequency 
of  the  generated  amplifying  resonance,  where  the  group  velocity  is  minimized.  Figure  SI 
shows  the  gain  exponent,  and  group  index  as  functions  of  frequency. 

BER  Calculation 

The  bit -error-rate  (BER)  is  determined  using  the  following  integrate-and-dump  matched- 
filter  procedure.  For  each  symbol  (0  and  1).  there  arc  50  test,  pulse  waveforms.  For  each 
test  waveform,  two  reference  waveforms  Ro[t)  and  R\(t)  are  generated.  Each  reference 
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waveform  Rj('t)  is  an  average  of  49  pulse  waveforms  of  symbol  type  j  and  dot's  not  include 
the  test  waveform. 

For  each  Rj(t),  the  normalization  integral 

Nj(r )  =  R2j(t)dt  (1) 

is  calculated,  where  ts  is  a  time  chosen  arbitrarily  on  the  leading  edge  of  the  pulse.  Next, 
the  normalization  constants 

**-^r*«*  (2) 

are  determined,  where  x(t)  is  the  test  pulse  waveform  and  ra  is  a  normalization  integration 
time  chosen  arbitrarily  on  the  leading  edge  of  the  Gaussian  pulse  before  the  point  where 
the  symbols  separate.  The  integral  D(r)  —  /o(r)  —  I\(r)  is  determined  for  each  test 
waveform,  where 

W)  =  ^  L  M  r+T  x(t)Rj(t)dt  ( j  =  0, 1).  (3) 

oij{Ta)^](r)  Jtm 

Each  integral  Ij(r)  is  a  measure  of  of  the  similarity  between  the  test  waveform  x(t )  and 
the  reference  waveform  Rj(t)  over  the  time  interval  from  t,  to  (tt  +  r). 

For  each  symbol  type,  the  50  test  pulse  waveforms  are  used  to  generate  a  probability 
density  of  D{r),  which  is  then  fit  to  a  Gaussian  distribution  normalized  to  an  area  of  1/2. 
The  overlapping  area  of  these  two  Gaussian  distributions  is  the  BER. 

Fig.  Si.  Gain  and  group  index,  (a)  The  dotted  line  shows  the  measured  gain  exponent 
(gL)  and  the  solid  line  shows  the  Lorentzian  least-squares  fit.  (b)  The  theoretically 
predicted  group  index  ng  corresponding  to  the  Lorentzian  gain  resonance  shown  in  (a). 
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Abstract:  We  describe  a  methodology  to  maximize  slow-light  pulse 

delay  subject  to  a  constraint  on  the  allowable  pulse  distortion.  We  show 
that  optimizing  over  a  larger  number  of  physical  variables  can  increase  the 
distortion-constrained  delay.  We  demonstrate  these  concepts  by  comparing 
the  optimum  slow-light  pulse  delay  achievable  using  a  single  Lorent/ian 
gain  line  with  that  achievable  using  a  pair  of  closely-spaced  gain  lines.  We 
predict  that  distortion  management  using  a  gain  doublet  can  provide  over  a 
factor  of  6  increase  in  slow-light  pulse  delay  as  compared  with  the  optimum 
single-line  delay.  Experimental  results  employing  Brillouin  gain  in  optical 
fiber  confirm  our  theoretical  predictions. 
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Over  the  last  decade,  there  has  been  great  progress  in  devising  new  methods  for  tailoring  the 
dispersion  of  opticaJ  materials,  such  as  electro  magnetically  induced  transparency  [I],  photonic 
crystals  [2],  and  nano-optic  resonators  [3-5],  This  work  has  been  motivated  by  the  need  for 
electronically  or  optically  controllable  pulse  delays  for  applications  such  as  optical  buffers, 
optical  memories,  and  signal  processing.  In  these  applications,  the  primary  requirements  for 
slow- light  pulse  delay  are  that  the  temporal  pulse  delay  />/  be  large  relative  to  the  pulse  width  r 
and  that  the  pulse  not  be  substantiaJly  distorted  (defined  below).  These  two  requirements  largely 
oppose  each  other,  with  large  delay  coming  at  the  cost  of  greater  distortion.  These  tradeoffs 
have  been  studied  in  simple  Lorentzian  systems  [6-8]  and  in  Doppler- broadened  media  [9] 
with  encouraging  results,  but  distortion  remains  a  major  limitation  to  the  usefulness  of  slow- 
light  pulse  delay. 

In  this  paper,  we  present  a  method  for  achieving  large  slow-light  pulse  delay  under  the  con¬ 
straint  that  the  distortion  does  not  exceed  a  particular  limit.  We  also  present  a  simple  distortion- 
managed  medium  based  on  dual  Lorentzian  gain  lines.  The  dispersive  properties  of  this  medium 
are  tailored  by  adjusting  the  spacing  between  two  adjacent  gain  lines.  We  predict  a  factor  of  6.25 
improvement  in  the  relative  slow-light  pulse  delay  for  the  two-gain-line  case  in  comparison  to 
the  standard  single-gain-line  case  for  the  same  distortion  constraint  in  both  cases.  We  demon¬ 
strate  dispersion-managed  pulse  delay  in  an  experiment  where  the  slow-light  effect  arises  from 
the  dispersion  associated  with  stimulated  Brillouin  scattering  (SBS)  resonances  in  a  pumped 
optical  fiber  [10-13]. 

For  pulses  propagat  ing  through  linear  optical  systems,  the  output  pulse  amplitude  A  (a),  r)  in 
the  frequency  domain  can  be  related  to  the  input  pulse  amplitude  A (co, 0)  by 

A{0>,i)  =/4(ft).0)e1*t“):J  (1) 

where  z  is  the  length  of  the  medium  and  k{(o)  is  the  complex  wavenumber  as  a  function  of 
frequency  (o .  A  pulse  propagates  undistorted  through  a  dispersive  material  when  k(co)  takes 
the  form 

k(a>)^ko-\-k\{m-($c)t  (2) 

where  co(  is  the  carrier  frequency  of  the  pulse  and  k\  is  real.  That  is,  the  pulse  shape  remains 
unchanged  and  the  only  effects  of  propagation  are  delay,  an  overall  phase  shift,  and  gain  or 
attenuation.  In  this  ideal  case,  the  pulse  delay  t#  is  equal  to  the  group  delay  tg  —  z(k\  -  l/c), 
where  c  is  the  speed  of  light  in  vacuum.  In  practice,  dispersive  media  do  not  satisfy  Eq.  (2) 
precisely,  but  have  higher-order  terms  in  the  Taylor  expansion  k(co)  0 kj((o  -  (ocy / /!, 

where  kt  d^k(co)/d(0^ .  These  higher-order  terms  can  lead  to  pulse  distortion,  or  a  change  in 
the  pulse  shape.  One  can  always  use  pulses  with  narrower  bandwidth,  and  thereby  improve  the 
linearity  of  &(co)  over  the  pulse  bandwidth,  but  this  results  in  longer  pulses  without  increasing 
the  delay.  It  is  sometimes  possible  to  simultaneously  increase  k\  and  decrease  the  pulse  band¬ 
width,  which  can  result  in  larger  relative  delays,  but  may  result  in  unacceptably  large  gain  or 
absorption  [6]. 


One  method  for  creating  delayed  pulses  with  minimal  distortion  is  to  reduce  the  effects  of  the 
higher-order  terms  in  the  Taylor  series  expansion  of  k(co).  The  effects  of  these  terms  have  been 
studied  in  simple  systems  [6-9, 13],  but  additional  degrees  of  freedom  can  be  used  to  eliminate 
or  balance  them.  We  propose  using  custom  slow-light  media  designed  to  minimize  the  effects 
of  distortion.  Creating  a  custom  dispersion  profile  k{co)  is  in  general  quite  difficult.  However, 
one  can  instead  create  a  custom  system  by  combining  multiple  simple  systems.  We  consider 
here  a  system  composed  of  two  Lorentzian  gain  lines;  a  gain  doublet,  as  shown  in  Fig.  1.  The 
wavenumber  of  such  a  system  is  given  by 


*(<»)  =  -«o  +  - 
c  z 


_ 1 _ + _ 1 _ ) 

co-(coo~8)+iY  (o  -  (fflb  +  5)  +  iy/  ' 


where  »q  is  the  background  refractive  index,  go/z  is  the  line-center  amplitude  gain  coefficient 
for  each  line,  y  is  the  linewidth,  and  25  is  the  separation  between  the  lines.  In  many  slow-light 
media,  the  dominant  source  of  distortion  Is  the  second-order  term  in  the  expansion  of  £(<o), 
k 2  [6].  For  the  gain  doublet,  we  have 


4^o  y2  3 5?  -  y2 

*  (P+f? 


As  can  be  seen  from  Eq.  (4),  &2  =  0  for  a  separation  of  S  ~  y/ >/X  thereby  eliminating  the 
lowest-order  distortion-causing  term.  For  a  doublet  in  this  configuration,  the  group  delay  is 
given  by 

=  =  +  (5) 

c  4  y 

which  should  be  compared  to  the  delay  for  a  single  line 


By  comparing  Eqs.  5  and  6,  it  is  seen  that  the  group  delays  for  these  two  systems  are  very  simi¬ 
lar.  While  this  result  is  interesting  and  illuminating,  it  is  of  limited  practical  value;  it  eliminates 
only  second-order  distortion,  ignoring  all  higher-order  terms. 

A  more  practical  approach  to  distortion  management  is  to  measure  distortion  as  the  deviation 
of  the  medium  from  the  ideal  one  described  in  Eq.  (2).  At  this  point,  it  becomes  convenient  to 
describe  the  medium  in  terms  of  its  transfer  function  //(ft)),  which  is  easily  related  to  k (ft)) 
by  //(ft))  exp(/A:(ft))z).  An  ideal  medium  has  a  transfer  function  with  constant  amplitude 

| //(ft));  =  Hq  and  a  phase  that  varies  linearly  with  frequency  lH[(o)  tp(Q ,  where  tp  is  the 
total  propagation  time  of  the  pulse.  Relating  these  to  Eq.  (2)  yields  Ho  -  exp!/(#o  k\  coo) z\  and 
tp  =  k\z.  Any  deviation  from  this  behavior  leads  to  pulse  distortion  [14]. 

There  are  many  ways  to  quantify  the  deviation  of  the  transfer  function  from  ideal.  One  very 
simple  way  is  to  calculate  the  infinity-norm — the  maximum  magnitude — of  the  amplitude  and 
phase  deviation  from  ideal.  This  leads  to  two  distortion  metrics,  one  for  the  amplitude  variation 
C Da)  and  one  for  the  phase  variation  ( Dp ).  The  amplitude  distortion  is  given  by 


(7) 


where  Hmm  and  //max  are  the  minimum  and  maximum  values  of  //(ft))  over  the  frequency 
range  (ft) q  —  A/,,  fth  4*  A^).  Similarly,  we  define  the  phase  distortion  as 


DP  =  —maxlUH(co)  -  (yj>+*>)|]J5!£, 


(8) 


frequency  (to  -  (oQ)/y 

Fig.  1 .  Gain  and  dispersion  for  a  single  Lorentzian  line  (solid)  and  a  doublet  with  separation 
5  -  yjs/i  (dashed).  Also  shown  are  the  two  constituent  lines  (dotted)  that  make  up  the 
doublet,  (a)  The  gain  exponent  has  a  broad  flat  top  (a  result  of  setting  £2  0).  although  it 

is  larger  for  the  same  value  of  the  individual  gain  coefficients,  (b)  The  dispersion  is  very 
similar  for  both  situations  near  the  center  of  the  lines,  but  extends  farther  for  the  gain- 
doublet  case. 


where  tp  and  $o  are  chosen  to  minimize  Dp.  This  provides  a  very  conservative  distortion  mea¬ 
sure  for  pulses  whose  power  falls  entirely  within  the  frequency  range  (Cty  -  A/,,  coo  -FA/;).  The 
calculation  of  Dp  also  serves  to  define  an  effective  propagation  time  tp  as  the  propagation  time 
of  a  pulse  through  the  ideal  medium  most  closely  approximated  by  the  real  medium.  A  delay 
defined  in  terms  of  this  propagation  time  t(i  =tp-  n^z/c  is  similar  in  concept  to  the  group  delay, 
except  that  it  is  based  on  the  material  dispersion  over  the  entire  bandwidth  of  interest  rather  than 
just  at  the  carrier  frequency.  In  fact,  the  pulse  delay  t(i  is  exactly  equal  to  the  group  delay  of  the 
best-fit  ideal  medium.  Because  this  delay  includes  dispersive  effects  over  the  entire  bandwidth 
of  interest,  it  provides  a  better  prediction  of  actual  pulse  delay  (peak  delay,  for  example)  than 
the  group  delay.  As  is  common  in  communication  systems  [14],  pulses  with  power  extending 
outside  this  range  can  still  be  used,  but  the  result  is  then  approximate.  These  distortion  met¬ 
rics  are  reasonable  in  the  case  where  the  signal  spectrum  is  approximately  limited  to  the  given 
bandwidth  window  but  is  otherwise  unknown.  If  more  is  known  about  the  signal  spectrum,  then 
more  complicated  distortion  measures  can  be  used  that  consider  a  fixed  spectrum  and  the  form 
of  the  material  deviation  from  ideal. 

Having  defined  these  distortion  metrics,  we  can  use  them  to  explore  the  gain-doublet  slow- 
light  media  discussed  above.  I  jet  us  first  consider  the  simple  Lorentzian  gain  line,  whose  trans- 


fer  function  can  be  written  as 


//it®) 


exp  ( izno ®)  xexp(j?1(ft))) 

/.  ©  /'y  \ 

exp  /zno—  +  #017 - r— —  • 

\  c  (a-03o)  +  iyj 


(9) 


Holding  y  fixed,  we  find  that  for  each  A$  there  is  an  optimum  value  of  got  that  provides  max¬ 
imum  delay  tc\  subject  to  gain  and  distortion  constraints.  We  have  chosen  to  use  an  amplitude 
gain  constraint  of  #qj  <  2.5  so  that  the  system  can  be  reasonably  implemented  experimentally. 
Larger  gain  is  often  accompanied  by  nonlinear  optical  effects. 

The  communication  system  of  interest  will  define  both  the  appropriate  form  of  the  relevant 
distortion  and  its  maximum  allowed  value.  In  conjunction  with  the  infinity  norm,  we  select  a 
conservative  value  of  maximum  allowed  distortion  Da  <  0.05  and  Dp  <  0.05.  This  value  is 
somewhat  arbitrary  and  results  in  pulse  distortion  that  is  just  noticeable  upon  visual  inspection. 

Figure  2(a)  shows  the  maximum  relative  delay  r^A/7  as  a  function  of  the  relative  bandwidth 
A n/y.  Figure  2(b)  show's  the  Lorentzian  gain  exponent  goi-  At  each  value  of  A/,,  the  gain  is 
chosen  to  maximize  the  delay  subject  to  the  gain  and  distortion  constraints  described  above. 
By  comparing  the  solid  lines  in  Figs.  2(a)  and  2(b),  we  find  that  the  delay  is  limited  by  the 
distortion  constraint  for  large  relative  bandwidths — A t?/7>  0.2:  driving  the  gain  higher  will 
produce  too  much  distortion.  For  small  bandwidths,  the  delay  is  limited  by  the  gain  constraint. 
The  maximum  relative  delay  is  achieved  at  some  intermediate  bandwidth — A^/y  «  0.2 — for 
which  both  distortion  and  gain  limits  are  met  simultaneously.  At  this  value,  we  find  a  maximum 
relative  delay  of  t^Af,  ^  0.5  under  the  constraint  that  neither  Dn  nor  Dp  exceed  0.05. 

We  now  consider  a  distortion-managed  system  constructed  using  two  nearby  Lorentzian  gain 
lines.  The  transfer  function  of  this  system  is 
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As  above,  we  can  optimize  over  the  tree  parameters — now  go2  and  8 — to  maximize  the  delay 
at  each  A/;  subject  to  the  same  constraints,  D(l  <  0.05,  Dp  <  0.05.  and  max  |#2(a>)]  £  2.5. 
The  dashed  line  in  Fig.  2(a)  shows  distortion-  and  gain-limited  relative  delay  as  a  function  of 
bandwidth  for  this  doublet  system. 

As  described  above,  the  second  order  distortion- — dominant  in  many  low-distortion  cases- 
can  be  eliminated  entirely  (ki  =  0)  by  setting  8  yj\J 3  — ♦  8/yzz  0.58.  One  might  expect 

this  result  to  provide  the  optimal  line  separation.  However,  as  shown  in  Fig.  2(c),  delay  can 
be  increased  at  large  bandwidths  by  using  a  larger  separation  because,  although  there  is  a  gain 
dip  at  the  carrier  frequency,  the  gain  excursion  over  the  full  bandwidth  can  be  minimized.  For 
small  bandwidths,  a  smaller  separation  is  desirable  because  the  delay  is  limited  by  gain,  not 
distortion,  and  so  it  is  more  useful  to  make  the  separation  small,  creating  larger  dispersion  for 
the  same  maximum  gain. 

There  are  three  distinct  frequency  regions  on  these  plots.  The  left-most  region  (A h/Y<  0.25) 
corresponds  to  gain-limited  configurations  where  the  second  line  provides  no  additional  benefit. 
The  right  most  region  (A b/Y>  0.8)  corresponds  to  the  distortion-limited  regime.  In  the  central 
region,  both  constraints  are  hit  simultaneously.  In  this  gain-doublet  medium,  we  observe  a  max¬ 
imum  relative  delay  of  0.9  at  a  bandwidth  of  Ah/y  —  0.8,  an  improvement  by  a  factor  of  6.25 
over  the  single  line  at  that  bandwidth  and  nearly  twice  the  best  single-line  relative  delay. 

To  test  our  distortion-management  concept,  we  propagate  optical  pulses  through  a  SBS  slow- 
light  medium  [10-13],  A  gain-doublet  with  an  adjustable  frequency  separation  28  is  realized 


Fig.  2.  Simulation  and  experimental  results  for  both  die  single  Lorentzian  line  (solid  lines 
for  simulation,  circles  for  experimental  results)  and  the  doublet  (dashed  lines  for  simula¬ 
tion,  squares  for  experimental  results),  (a)  Relative  delay,  (b)  Lorentzian  line-center  ampli¬ 
tude  gain  coefficients  #o  1  and  £q2.  Also  shown  is  the  center- frequency  gain  coefficient  for 
the  doublet  g2(^)  (dot-dashed),  (c)  Line  separation  for  the  doublet.  The  horizontal  line 
indicates  the  value  of  d/y  that  leads  to  *2  —  0 


TL1  IS1  FPC1  MZM1  Cl 


Fig.  3.  Experiment  setup  based  on  a  fiber  BriUouin  amplifier.  TL 1 .  TL2:  tunable  lasers:  IS  I . 
IS2:  isolators:  FPC 1,  FPC2.  FPC3:  fiber  polarization  controllers;  MZM1.  MZM2:  Mach- 
Zehnder  modulators;  FG1,  FG2:  function  generators;  EDFA:  Erbium  doped  fiber  amplifier: 
C 1 ,  C2:  circulators;  SMF-28e:  500-m-long  SMF-28e  fiber  (the  SBS  amplifier);  PM:  power 
meter. 


by  pumping  a  standard  room -temperature  telecommunications  optical  fiber  with  a  bichromatic 
laser  field,  which  is  obtained  via  cairier-frequency -suppress ion  modulation  of  a  continuous- 
wave  (CW)  single- frequency  laser  field  [15].  The  doublet  spacing  25  is  adjusted  by  changing 
the  frequency  of  the  voltage  driving  the  modulator  and  the  gain  coefficient  (and  hence  the  slow- 
light  delay)  is  ad  justed  by  changing  the  pump  power.  Either  of  these  adjustments  can  be  made 
very  quickly,  limited  by  the  transit  time  of  the  pump  field  through  the  optical  fiber. 

The  transfer  function  //2(G))  (Eq.  10)  can  be  realized  by  setting  the  carrier  frequency  of  the 
pulses  close  to  the  SBS  amplifying  resonances  (the  .so-called  Stokes  resonances)  and  keeping 
the  pulse  intensity  low  enough  so  that  the  response  is  linear.  For  comparison,  we  also  propagate 
pulses  through  a  standard  slow-light  medium  with  a  single  gain  resonance,  corresponding  to 
H\(co)  (Eq.  9)  and  realized  using  a  single-frequency  pump  field  [10-13]. 

The  experimental  setup  is  shown  schematically  in  Fig.  3.  Two  1550-nm  narrow-linewidth 
tunable  lasers  (TLI  and  TL2)  are  used  to  produce  Stokes  signal  pulses  and  the  pump  beams, 
respectively.  A  bichromatic  pump  beam  is  produced  by  passing  the  single -frequency  CW  laser 
beam  generated  by  TL2  (angular  frequency  (op)  through  modulator  MZM2,  which  is  driven  by 
a  sinusoidal  voltage  of  angular  frequency  5,  The  bias  voltage  applied  to  MZM2  is  chosen  to 
suppress  the  carrier  frequency,  resulting  in  a  beam  with  frequencies  (0)p  T  5).  This  beam  is  am¬ 
plified  by  an  erbium-doped  fiber  amplifier  (EDFA)  and  routed  via  circulator  C2  to  pump  a  500- 
m-long  SMF-28g  fiber  (the  slow-light  medium).  For  this  fiber,  we  measure  an  SBS  linewidth 
(full-width  at  half-maximum)  of  y//r  ^  35  MHz.  After  passing  through  the  slow-light  medium, 
the  beam  is  routed  out  of  the  system  via  circulator  Cl  and  monitored  by  a  power  meter. 

The  signal  pulses  are  produced  from  a  laser  beam  generated  by  TLI.  The  beam  passes 
through  an  isolator,  a  fiber  polarization  controller,  and  modulator  MZM1,  producing  pulses  with 
carrier  frequency  coc.  They  are  approximately  Gaussian -shaped  with  an  intensity  envelope  of 
the  form  l{t)  =  /oexp  - (/ / r)2],  where  the  spectrum  S(eo)  expi  r2(a>  -  a)f)2]  has  a  1  je  half 
width  of  Aq  1  /t.  The  pulses  enter  the  500-m-long  SMF-28e  fiber  via  circulator  C I  and  coun- 
terpropagate  with  respect  to  the  bichromatic  pump  beam  in  the  fiber.  The  sJow-light-delayed 
and  amplified  pulses  are  routed  out  of  the  system  via  C2,  sensed  by  a  fast  photodetector,  and 
displayed  on  a  digital  oscilloscope.  The  frequency  difference  between  TLI  and  TL2  is  tuned 
so  that  the  pulse  carrier  frequency  is  set  precisely  to  the  center  of  the  SBS  amplifying  reso¬ 
nances,  (Oc  —  cqq  where  a*)  =  cop  —  Q/j,  and  0.^j2n  1 2.5  GHz  is  the  Brillouin  frequency  shift 


for  rhe  SMF-28e  fiber.  Fiber  polarization  controllers  FPC1  and  FPC2  are  used  to  maximize  the 
transmissions  through  the  Mach-Zehnder  modulators,  and  FPC3  is  used  to  maximize  the  SBS 
slow-light  delay  experienced  by  the  Stokes  pulses. 

In  the  experiments  with  the  distortion-managed  slow-light  medium,  we  select  a  pulse  band¬ 
width  A i,  (using  the  assumption  that  =  Ao),  set  the  resonance  half- separation  S  according 
to  Fig.  2(c),  set  the  central  gain  according  to  Fig.  2(b)  by  adjusting  the  pump  power, 

and  measure  the  delay  of  the  peak  of  the  pulses  induced  by  the  slow-light  medium  (which  is 
approximately  equal  to  / <j).  From  Fig.  2(a),  it  is  seen  that  the  agreement  with  the  theoretical 
predictions  is  very  good.  Also  shown  in  Fig.  2(a)  are  our  observations  for  a  single-gain  line  (the 
non-distortion-managed  slow-light  medium),  where  the  agreement  between  the  observations 
and  the  predictions  is  excellent. 

The  data  shown  in  Fig.  2(a)  demonstrate  that  the  distortion-managed  slow-light  medium 
vastly  outperforms  the  standard  slow-light  medium  under  conditions  of  constant  distortion. 
Around  the  optimum  normalized  bandwidth  of  ~0.84  (r  ^  1 1  ns),  we  measure  an  improve¬ 
ment  in  relative  pulse  delay  of  a  factor  of  8.1  ±  2.6.  Thus,  the  distortion-managed  approach 
substantially  increases  the  usable  bandwidth  of  a  slow-light  medium  that  is  much  easier  to 
implement  in  comparison  to  previously  suggested  methods  [16, 17]. 

In  conclusion,  we  have  shown  that,  although  distortion  limits  the  slow-light  delay  achievable 
in  simple  systems  [<S-8J,  it  is  possible  to  achieve  better  distortion-  and  gain-limited  delay  using 
custom  composite  media.  Dramatic  improvement  can  be  achieved  over  a  single  Lorentzian  line 
with  a  system  as  simple  as  a  Lorentzian  doublet,  suggesting  that  even  greater  improvement  may 
be  achievable  using  more  flexible  media. 

We  gratefully  acknowledge  the  financial  support  of  the  DARPA  DSO  Slow-Light  Program. 
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The  causality  of  waves  propagating  in  a  double-negative  (DNG)  meta  material  (e,.<0  and  yu,,  <0)  has  been 
investigated  both  analytically  and  numerically.  By  considering  the  one-dimensional  electromagnetic  problem 
of  a  pulsed  current  sheet  radiating  into  a  DNG  medium,  it  is  shown  that  causality  is  maintained  in  the  presence 
of  a  negative  index  of  refraction  only  if  the  DNG  medium  is  dispersive.  A  Drude  model  DNG  medium  is  used 
in  this  study.  Spectrograms  of  the  wave  phenomena  in  the  dispersive  DNG  medium  show  that  the  higher 
frequency  components,  which  create  the  leading  edge  of  the  electromagnetic  signals  and  sec  a  double  positive 
(DPS)  medium  (er>0  and  pLr> 0),  arrive  causally  before  the  negative  index  effects  germinate  completely. 
Comparisons  with  approximate  analytical  results  demonstrate  the  presence  of  the  negative  index  of  refraction 
properties  in  the  continuous  wave  portion  of  the  signals.  This  dynamic  pulse  reshaping  between  the  positive 
and  negative  index  of  refraction  wave  components  causes  an  apparent  delay  in  the  realization  of  the  negative 
index  of  refraction  properties.  Pulse  broadening  of  the  signal  tails  is  associated  with  both  dispersion  and  a 
larger  negative  index  of  refraction  seen  by  the  associated  wave  components. 
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I.  INTRODUCTION 

The  index  of  refraction  of  a  double- negative  (DNG) 
metamaterial,  i.e.,  a  material  with  both  negative  permittivity 
and  negative  permeability  (see,  for  instance,  [1-4]),  has  been 
shown  to  be  negative  (see,  for  instance,  [5-7]).  There  are 
now  several  theoretical  and  experimental  studies  that  have 
been  reported  confirming  this  negative  index  of  refraction 
(NIR)  property  and  applications  derived  from  it,  such  as 
phase  compensation  and  electrically  small  resonators  [8], 
negative  angles  of  refraction  [8-13],  enhanced  focusing 
[14-17],  backward  wave  antennas  [18],  Cerenkov  radiation 
[19],  photon  tunneling  [20,21],  and  enhanced  electrically 
small  antennas  [22].  These  studies  rely  heavily  on  the  con¬ 
cept  that  a  continuous  wave  (CW)  excitation  of  a  DNG  me¬ 
dium  will  lead  to  a  NIR  and,  hence,  to  negative  phase  terms. 

If  one  considers  carefully  the  ramifications  of  a  homoge¬ 
neous,  nondispersive  DNG  medium  and  the  resulting  NIR, 
one  immediately  encounters  a  paradox  in  the  time  domain.  A 
source  in  such  a  DNG  medium  will  generate  a  signal  that 
propagates  noneausally.  On  the  other  hand,  many  of  the  NIR 
studies  have  been  based  on  numerical  simulations  in  the  time 
domain  of  dispersive  DNG  meta  mate  rials.  In  fact,  as  shown 
in  [5],  realistic  DNG  metamaterials  must  be  dispersive.  How 
then  does  one  reconcile  the  CW  NIR  properties  associated 
with  a  nondispersive  DNG  medium  with  causa)  pulse  propa¬ 
gation  in  a  dispersive  DNG  medium?  Can  one  rely  on  de¬ 
signs  of  NIR  applications  that  rely  on  the  CW  behavior  of 
the  corresponding  nondispersive  DNG  medium? 

In  this  paper,  the  issue  of  causality  in  DNG  metamaterials 
is  considered  by  studying  the  one-dimensional  electromag¬ 
netic  problem  of  a  pulsed  source  in  a  DNG  medium.  A  lossy, 
dispersive  Drude  model  is  used  to  create  the  DNG  medium. 
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A  finite  difference  time  domain  (FDTD)  numerical  solution 
of  the  problem  is  obtained.  The  source  is  driven  by  an  exci¬ 
tation  pulse  that  is  based  on  a  sinusoid  and  a  smooth  enve¬ 
lope.  It  is  demonstrated  that  both  the  expected  causal  propa¬ 
gation  and  the  NIR  effects  are  obtained.  Other  recent 
considerations  of  causality  (e.g..  [13,23])  have,  in  effect, 
been  limited  to  narrow  bandwidth  considerations.  The  corre¬ 
sponding  nondispersive  DNG  medium  problem  is  then  con¬ 
sidered.  It  is  demonstrated  analytically  that  only  a  noncausal 
exact  solution  exists.  The  noncausal  analytical  DNG  .solution 
and  the  causal  FDTD  DNG  results  are  then  reconciled  by  the 
presence  of  dispersion.  It  is  shown  that  the  leading  edge  of 
the  pulse,  which  is  associated  with  high  frequencies,  propa¬ 
gates  causally  in  a  double-positive  (DPS)  medium  (i.e.,  a 
normal  medium  that  has  both  positive  permittivity  and  posi¬ 
tive  permeability),  while  the  CW  portion  of  the  pulse  senses 
the  DNG  medium  and  exhibits  the  NIR  effects.  The  slower 
speed  components  appear  to  accumulate  in  the  trailing  edge. 
Both  the  leading  and  trailing  edges  of  the  pulse  are  shown  to 
be  strongly  affected  by  the  presence  of  dispersion.  Using  the 
nondispersive  results  to  guide  the  construction  of  approxi¬ 
mate  solutions  to  the  dispersive  DNG  medium  problem,  it  is 
shown  that  the  NIR  effects  obtained  for  the  dispersive  case 
coincide  with  those  predicted  by  the  nondispersive  case.  It  is 
also  shown  that  observed  time  delays  in  the  realization  of  the 
NIR  properties  are  directly  associated  with  the  transition 
from  the  causal  DPS  response  to  the  DNG  one.  These  obser¬ 
vations  reinforce  the  fact  that  DNG  metamaterials  and  the 
associated  NIR  results  are  physical  and  that  their  experimen¬ 
tal  confirmation  is  consistent  with  the  underlying  wave 
propagation  physics. 

11.  1D-FDTD  SIMULATOR 

A  one-dimensional  FDTD  simulation  environment  was 
utilized  to  study  the  characteristics  of  signals  produced  by  a 
pulsed  source  in  a  DNG  medium.  This  one-dimensional  en~ 
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vironrnent  was  convenient  and  incorporated  all  of  the  neces¬ 
sary  physics.  The  held  components  were  assumed  to  be  Ex 
and  Hv  with  the  direction  of  propagation  taken  along  the  z 
axis.  An  electric  current  plane-wave  source  was  located  in 
the  plane  z  =  0,  at  the  center  of  a  very  large  DNG  slab.  The 
DNG  slab  was  surrounded  by  free  space  and  was  centered  in 
the  FDTD  domain.  This  slab  geometry'  was  used  in  order  that 
standard,  exact  absorbing  boundary  conditions  could  be  ap¬ 
plied  at  the  edges  of  the  FDTD  simulation  space.  The  prob¬ 
lem  geometry  is  shown  in  Fig.  1. 

As  in  [6],  lossy  Drude  polarization  and  magnetization 
models  were  used  to  simulate  the  DNG  medium.  In  the  fre¬ 
quency  domain,  assuming  an  exp(  —  icot)  time  dependence, 
this  means  the  permittivity  and  permeability  were  described 
as 


e(w)  =  e0  1  - 


p* 


w(c 0  +  iTe)j  * 


mM  =  mo!  * 


to 


pm 


io(u)+  lTm)  I 


(I) 


The  corresponding  time-domain  equations  for  the  polariza¬ 
tion,  Px,  and  the  normalized  magnetization,  Mn^M y 
fields  were 


d;Px+r ,<?,/>,= e0topeEx , 

r>]M ny  +  ^ mdtM ny = •  (2) 

The  normalized  magnetization  was  introduced  to  make  the 
electric  and  magnetic  held  equations  completely  symmetric. 
By  introducing  the  induced  electric  and  magnetic  currents 


Jx-Wx. 


Ky  =  d,M 


nv  » 


(3) 


*  » 


Current  sheet 


FIG.  1.  (Color  online)  Configuration  of  the  one-dimensional 
electromagnetic  current  sheet  problem. 

1 

d(Ex= - (^//v  +  iv), 

eo 

1 

djHx— - dzEx .  (5) 

Mo 

Equation  sets  (4)  and  (5)  were  solved  self-consistent ly 
and  numerically  with  the  FDTD  approach  [24.25],  i.e.,  these 
equations  were  discretized  with  a  standard  leap-frog  in  time, 
staggered  grid  approach.  The  electric  field  component  was 
taken  at  the  cell  edges  for  integer  time  steps;  the  magnetic 
field  components  were  taken  at  the  cell  centers  for  half- 
integer  time  steps.  The  electric  and  magnetic  currents  were 
located  in  space  at  the  cell  centers,  and  their  time  assign¬ 
ments  were  opposite  to  the  corresponding  electric  and  mag¬ 
netic  field  components,  i.e.,  the  magnetic  current  components 
were  sampled  at  integer  time  steps  and  the  electric  current 
component  was  sampled  at  half-integer  time  steps.  This  al¬ 
lowed  a  FDTD  stencil  that  properly  simulated  matched  me¬ 
dium  conditions.  The  time  step  was  set  at  0.95  of  the  Courant 
value,  i.e.,  A/=0.95Az/c. 

In  all  cases,  the  center  frequency  of  interest  was  chosen  to 
be  /0s£  <u0/27r=30  GHz,  corresponding  to  a  fnee-space 
wavelength  \o=1.0cm.  Only  matched  media  at  f0  were 
considered,  i.e., 


and  the  source  current,  the  field  and  current  equations 
used  to  model  the  waves  generated  by  a  pulsed  source  in  a 
DNG  medium  became 


Z(«q)  - 


//*(<■>()) 
V  e(<i>0) 


(6) 


l 

d(Ex— - ( d:Hy  +  Jx+  Js)t 

so 

<Vjt  +  ^  x~  ^0  * 

l 

d'Hv= - (dtEx  +  Ky), 

Mo 

d,Ky+  r mKy  =  /j.0(o2pmHy .  (4) 

A  lossless,  nondispersive  DPS  medium  was  used  as  a  com¬ 
parison  case.  In  such  a  DPS  medium,  the  corresponding  field 
and  current  equation  set  was  simply 


Matching  was  obtained  by  setting  the  parameters  for  the 
electric  and  magnetic  Drude  models  to  be  identical,  i.e., 
(ope=  a>pitl=  a)p  and  r,=  rm  =  r.  In  all  cases,  only  low  loss 
values  were  considered  by  setti  ng  F  =  1 0 + 8  =  5 . 3 1 
X  10~4o>0.  This  value  was  selected  to  connect  these  results 
to  those  presented  in  [6]  and  [12].  A  DNG  medium  matched 
to  free  space  at/0=30  GHz  was  considered,  i.e., 

/eOo)  /m("o)  - 7  I — ; - 7 

BW-V';  "  v— — =  VM«o)v/M«o)  =  1. 


This  required  &>/,=  27rv2/0=2.66573X  10"  and,  hence,  F 
=  3.75X  10"4w/). 

The  source  was  given  by  the  expression 
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J(z)=Js{z)X=~S(z)S(t)x.  (8) 

It  was  modeled  in  the  FDTD  simulation  space  by  introducing 
the  equivalent  source  7Ji=|X//yv,  i.e.,  Hv$(mcc=S(t)  at 
the  source  plane.  The  input  time  signal  was  causal  and  was 
defined  by 

=  (9) 


where  the  multiple  cycle  m-n-m  envelope  is  given  by  the 
expression 


r  0 

for  /<0 

for  0^/< m Tp 

g(l] 

1 

to  r  rn  Tp**t<(  m  +  n )  Tp 

£  off ( * ) 

for  (m  +  n )  Tp ^  / <  (m  +  n  T  m )  Tp 

,0 

for  ( m  +  n  +  m )  Tp  ^  t , 

(10) 

where 

the 

period 

Tp=  l//().  With  the  terms  .vun(/) 

=  //(/// Tp)  and  =  —  the  continu¬ 

ous,  two  derivative  smooth  functions  gon  and  g0II  can  be 
written  as 

*«,(*)=  104,(0- 154,(0+64,(0. 

««n<0=  l  “[  104n(0-  15xJo(0+64fl(0]-  (1  0 

The  function  gm(t)  goes  smoothly  from  0  to  1  in  m  peiiods; 
the  function  gaff(t)  goes  smoothly  from  1  to  0  in  m  periods. 
Tile  function/!/)  thus  turns  on  smoothly  in  rn  periods,  main¬ 
tains  a  constant  amplitude  for  n  periods,  and  turns  off 
smoothly  in  m  periods.  The  more  cycles  in  the  center  por¬ 
tion,  i.e.,  the  larger  n  is,  the  narrower  the  bandwidth  and  the 
closer  the  input  pulse  is  to  a  CW  signal.  The  shorter  the 
turn-on  and  turn-off  sections  are,  i.e.,  the  smaller  m  is,  the 
broader  the  bandwidth  is  and,  in  particular,  the  more  high 
frequency  content  there  is  in  the  pulse. 

The  simulation  space  was  discretized  into  cells  with  a 
length  A;:  =  X.o/100=  100  /tm,  The  FDTD  predicted  electric 
field  values  at  c=  100 Ac  and  z~  120 Ac  for  a  5-10-5  input 
pulse  in  the  dispersive  DNG  medium  are  shown  in  Fig.  2(a). 
An  expanded  view  of  the  leading  edge  of  these  pulses  is 
shown  in  Fig.  2(b).  The  corresponding  FDTD  predicted  elec¬ 
tric  field  values  for  a  lossless  nondispersive  DPS  medium 
with  n—  +  I  are  shown  in  Figs.  3(a)  and  3(b).  Comparing 
these  sets  of  figures,  one  can  clearly  see  that,  as  in  the  DPS 
case,  the  fronts  of  the  signals  in  the  DNG  medium  arrive 
causally  at  the  observation  points.  In  the  DPS  case,  the  peaks 
of  the  central,  CW  portion  of  the  pulse  observed  at  the  more 
distant  point  2  occur  later  in  time  than  they  do  at  point  7.  In 
contrast,  one  can  see  that  in  the  DNG  case  the  NIR  effects 
(e.g.,  negative  phase  and  phase  velocity)  evolve  as  the  CW 
portion  of  the  signal  is  realized,  e.g.,  the  peaks  of  the  CW 
portion  of  the  pulse  observed  at  the  more  distant  point  2 
occur  earlier  in  time  than  they  do  at  point  /.  Dispersive 
effects  are  also  apparent  in  the  DNG  results  as  the  observa¬ 


FIG.  2.  (Color  online)  FDTD  predicted  electric  field  values  in 
the  DNG  medium  for  the  5-10-5  input  pulse  (a)  measured  at  the 
points  c  =  1 00 A z  and  z  —  1 20A z  and  (b)  zoomed  to  the  front  of 
those  pulses,  where  A z=  100  /ttm. 

tion  points  move  farther  away  from  the  source.  These  results 
are  consistent  with  those  reported  in  [6]  for  a  DNG  slab. 

III.  ANALYTICAL  CONSIDERATIONS 

Consider  now  the  plane- wave  source  in  a  hypothetical 
lossless,  nondispersive  DNG  medium  with  er=  pLr=  —  1 .  The 
Maxwell  equations  become 

i 

d,Hy=y-\dzEx.  (12) 

I I 

First,  let  the  excitation  signal  be  given  by  the  expression 

S(t)  =  [HU)-mt~nisin(<ont),  (13) 
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Point  l 

- 

sgn(;)  = 

•  Point  g 

-  1  for  z<0 
+  1  for  1 3*0. 


(17) 


FIG.  3.  (Color  online)  FDTD  predicted  electric  field  values  in 
the  DPS  medium  for  the  5- 10-5  input  pulse  (a)  measured  at  the 
points  z  =  lOQAz  and  z  —  I20Az  and  (b)  zoomed  to  the  front  of 
those  pulses,  where  Az=  100  /u,m. 


where  T={m~\~n  F  ?n)Tn  and  the  Heaviside  function 


H{ty 


0  for  f<0 
1  for  f^0. 


Introducing  the  variable 

'o=|e 


(14) 


(15) 


We  note  for  completeness  that  <?,|z|  =  sgn(z)  and  d ,  sgn(z)= 

+2  %z). 

The  corresponding  DPS  medium  results  with  er~  fxr~ 
4- 1  are  described  by  the  Maxwell  equations 


^Er=-7^|(<LWN-fyy), 


'•//'  '  \pt\'1:EX' 


(18) 


and  tlieir  solutions 


*  \n\'n, 

EAzrO  =  2  {HU ~  to) ~H[r~( t0  +  T) ]} 

Xsin  [«0(f-/0)], 


«y(z.0=  jsgn(z){/#(/- tn)-H[t-Vo+T)]} 


X  sin[ft»o(/~ /q)]- 


(19) 


Comparing  the  two  results,  Eqs.  (16)  and  (19),  one  can 
clearly  see  that  the  lossless,  nondispersive  DNG  medium  so¬ 
lution  is  noncausal.  This  observation  is  in  complete  agree¬ 
ment  with  those  made  in  [5]. 

Similarly,  with  the  m-n-m  excitation  pulse  given  by  Eq. 
(9),  the  DNG  medium  solution  is 


I A 


t/2 


Ex(zj)  —  ^  tf(f+*o) sin[ wo( t “F ] » 


Hy(z,t)=*  rsgnU>g(f  +  /0)sin[w0(;  +  /0)J,  (20) 


and  the  DPS  medium  solution  is 


i  M 


1/2 


Ex(z,t)=  ~  T-TT7r«('-fo)s'n[«>o(,-,o)l. 


it  is  readily  shown  that  the  solution  to  Eqs.  (12)  and  (13)  is 

1  \n\m 

Ex( Z,t )  =  2  J^pTT  I  H( t  +  to)  -  H{ t  + 10- T) ]sin[  io0(t + 10 ) ] , 

Hy(z,t)=^sgn(z)[H(t+t0)-H(t  +  t0-T)] 

Xsin[w0(f+f0)J,  (16) 

where 


1 

Hy(z,t)=  ^sgn(z)£(r-f0)sin[wo('-'o)J-  (21) 

The  DNG  medium  result  is  again  noncausal  while  the  DPS 
medium  result  is  causal.  The  FDTD  results  shown  in  Fig.  3 
recover  the  DPS  solution  (21)  to  within  0.1%. 

Led  by  the  observed  combination  of  causal  envelopes  and 
NIR  effects  of  the  central  portions  of  the  DNG  FDTD  re¬ 
sults,  we  consider  the  following  approximate  analytical  solu¬ 
tions  to  the  fields  produced  by  the  abrupt  and  smoothed  win¬ 
dowed  sine  waveforms,  respectively: 
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FIG.  4.  (Color  online)  Approximate  electric  field  values  in  the 
DNG  medium  for  the  5-10-5  input  pulse  and  A z  =  100  /xm  mea¬ 
sured  at  the  point  z=  lOOAz. 

I  \n\m 

Xsin[w0(/+  /(,)], 

Xsin[«0(f+fo)].  (22) 

and 

l  \/x\m 

Ex(z,t)  =  j  -  fu)sin[o)0(f +/(,)], 

I 

Hy(z,t)=  -sgn(z)g(/-r0)sin[£a0(r+/0)]-  (23) 

These  solutions  would  allow  for  forward  propagation  in  time 
while  maintaining  the  solution's  negative  CW  phase  charac¬ 
teristics.  Inserting  the  proposed  solution  (22)  into  Eqs,  (12), 
one  finds  that  it  does  not  exactly  solve  them.  The  defects  in 
this  proposed  solution  occur  as  ^-function  contributions  at  t 
=  r0  and  t  =  t0+T.  Similarly,  inserting  the  proposed  solution 
(23)  into  Eqs.  (12),  one  finds  that  it  also  is  not  an  exact 
solution  and  its  failings  as  an  exact  solution  occur  as  defects 
proportional  to  the  derivatives  of  the  turn-on  and  turn-off 
portions  of  the  envelope  function  g.  Note,  however,  that  die 
CW  portions  of  both  of  these  proposed  solutions  do  satisfy 
the  Maxwell  equations.  Consequently,  since  the  defects  are 
highly  localized,  we  will  take  Eqs.  (22)  and  (23)  as  approxi¬ 
mate  solutions  for  the  DNG  case. 

The  electric  field  components  of  the  approximate  analyti¬ 
cal  solutions  (22)  and  (23)  for  the  5-10-5  excitation  pulse  are 
compared  at  the  point  2=  lOOAz  in  Fig.  4.  Approximation  1 
denotes  the  electric  field  component  given  in  Eqs.  (22);  ap¬ 
proximation  2  denotes  the  electric  field  component  given  in 
Eqs,  (23).  Complete  agreement  of  the  CW  portions  of  these 
results  is  observed,  as  expected.  Their  differences  oc- 
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FIG.  5.  (Color  online)  Comparison  of  the  FDTD  predicted  and 
approximate  analytical  electric  field  values  in  the  DNG  medium  for 
the  5-10-5  input  pulse  measured  at  the  points  (a)  z  =  3()A.:  and  (b) 
z  =  lOOAc,  where  Az  =  100  /txm. 

cur  simply  in  how  each  of  them  is  turned  on  and  off.  We  nole 
that  the  approximate  solution  (22)  contains  extremely  high 
frequency  content  because  of  its  abrupt  tum-on  and  turn-off. 
The  n-m-n  form  is  more  realistic  and,  as  Fig.  4  shows,  it 
contains  the  essential  characteristics  of  this  limiting  case. 

The  electric  field  component  of  the  approximate  solution 
given  by  Eqs.  (23)  for  the  5-10-5  excitation  pulse  is  com¬ 
pared  to  the  corresponding  FDTD  results  for  the  dispersive 
DNG  medium  at  the  points  2=  30 Ac  and  at  z=l0()Az  in 
Figs.  5(a)  and  5(b),  respectively.  The  corresponding  results 
for  a  2-16-2  excitation  pulse  were  also  obtained.  Because  the 
turn-on  and  tum-off  portions  of  the  2-16-2  pulse  are  shorter, 
the  time  discretization  was  halved  by  decreasing  the  spatial 
discretization  by  a  factor  of  2  to  Az/  =  5.0X10  5  m.  The 
electric  field  component  of  the  approximate  solution  given 
by  Eqs.  (23)  for  the  2-16-2  excitation  pulse  is  compared  to 
the  FDTD  result  for  the  DNG  medium  at  the  points,  z 
=  60A z'  and  z  =  200Azf ,  in  Figs.  6(a)  and  6(b),  respectively. 
These  observation  points  are  located  at  the  same  distances 
away  from  the  source,  z  —  0.003  m  and  z  =  0.01  m,  as  they 
are  for  the  5-10-5  excitation  pulse  case. 

From  both  sets  of  results,  one  can  clearly  discern  the 
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FIG.  6.  (Color  online)  Comparison  of  the  FDTD  predicted  and 
approximate  analytical  electric  held  values  in  the  DNG  medium  for 
the  2-16-2  input  pulse  measured  at  the  points  (a)  z=60Az'  and  (b) 
z  =  200A<-',  where  Ac' =  50  fin\. 

causal  nature  of  the  entire  FDTD  solution.  One  can  aJso  see 
that  the  NIR  effects  in  the  CW  portions  of  the  FDTD  and  the 
approximate  results  agree  quite  well.  Thus,  the  NIR  effects 
associated  with  a  CW  analysis  are  indeed  recovered  by  the 
lossy,  dispersive  FDTD  results.  Clearly,  if  the  loss  were 
larger,  the  clarity  of  these  CW"  components  will  be  lost 
among  the  dispersive  effects.  This  occurs  even  for  the 
present  cases  as  the  observation  points  move  much  farther 
away  from  the  source  plane. 

One  can  also  immediately  see  from  Figs.  5  and  6  that  the 
main  differences  between  the  FDTD  and  the  approximate 
results  occur  principally  in  the  tum-on  and  turn-off  portions 
of  the  envelope,  as  expected  from  the  approximate  solution 
analysis.  The  presence  of  dispersion  is  clearly  seen  in  these 
regions.  The  faster  turn-on  envelope  of  the  2-16-2  pulse 
shows  less  dispersive  effects  than  does  the  5-10-5  excitation 
pulse.  This  might  be  expected  since  the  turn-on  portion  of 
the  envelope  of  the  2-16-2  pulse  has  higher  frequency  con¬ 
tent  than  does  the  corresponding  portions  of  the  5-10-5  en¬ 
velope.  Consequently,  according  to  the  Drude  model,  it  sees 
a  medium  more  like  free  space  and  less  dispersion.  On  the 
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other  hand,  there  are  larger  differences  in  the  turn-off  por¬ 
tions  of  the  2-16-2  excitation  pulse  results  than  there  are  in 
those  produced  by  the  5-10-5  excitation  pulse.  This  trailing 
edge  has  more  of  its  own  higher  and  lower  frequency  com¬ 
ponents  that  are  subject  to  dispersion;  and  it  accumulates 
more  of  the  slower  moving,  lower  frequency  components 
generated  in  the  earlier  portions  of  the  pulse.  Moreover,  com¬ 
paring  the  set  of  results  at  the  point  closer  to  the  source  with 
die  one  farther  away  from  it,  one  observes  that  there  are  time 
delays  between  when  the  CW  NIR  effects  are  completely 
present  and  when  they  lose  their  dominance  to  the  dispersive 
effects  and  that  these  delays  increase  as  the  observation  point 
moves  farther  away  from  the  source. 


IV.  SPECTROGRAM  CONSIDERATIONS 


Flow  is  it  then  that  the  dispersive  effects  reconcile  the 
NIR  and  causal  behaviors?  To  address  this  issue,  we  have 
studied  spectrograms  [26,27]  of  the  source-generated  pulse 
as  it  arrives  at  various  observation  points  in  the  DNG  me¬ 
dium.  We  recall  that  the  magnitude  of  die  phase  and  group 
speeds  of  the  wave  in  the  low  loss  limit  are  given,  respec¬ 
tively,  by  the  expressions 

i  _  c  _  c  c°l 

W"’1  1 1  +  Re(  v'  ' 

(24) 


Vg((*))  = 


4-Re(^)  +  o)r9J[Re(^)]  or -Ecu', 


7<\  (25) 


where  the  real  part  of  the  susceptibility  Re(^)-Re 
[  —  «-F i.T )]***  — cajja?'.  We  also  recall  that  the  point  in 

the  frequency  spectrum  of  the  Drude  model  that  isolates  the 
DNG  and  the  DPS  regions  occurs  when  i.e., 

when  (o=cor  .  hence  when  /=  v2/()=  42.42  GHz.  Well 
above  this  frequency,  the  permittivity  and  permeability  are 
those  of  free  space.  Thus,  the  highest  frequency  components 
propagate  causally  in  a  medium  that  looks  like  free  space.  In 
particular,  these  components  are  the  ones  mainly  responsible 
for  the  tum-on  and  turn-off  portions  of  the  envelope.  The 
lower  frequencies  propagate  with  speeds  slower  than  their 
value  at  the  CW  portion  of  the  pulse.  w(),  i.e.,  |l?p((u0)| 

To  illustrate  these  effects,  the  spectrograms  of  the  1-1-1. 
2-16-2,  and  5-10-5  generated  electric  field  pulses  at  z 
=  0.01  m  in  the  dispersive  DNG  medium  were  obtained. 
They  are  plotted  in  Figs.  7-9,  respectively.  The  spectrograms 
were  generated  by  taking  Blackman- windowed  Fourier 
transforms  at  every  fourth  sample  in  the  signal  data  1 26.27J. 
The  number  of  signal  samples  analyzed  by  each  transform 
was  detenu ined  such  that,  after  the  data  selection  was  win¬ 
dowed,  the  selection  would  represent  at  least  a  few  and  at 
most  several  wavelengths,  depending  on  its  frequency  com¬ 
ponents.  As  a  result,  512  signal  samples  were  used  in  the 
transforms  of  the  5-10-5  pulse  results,  and  1024  signal 
samples  were  used  in  transforms  of  the  1-1-1  and  2-16-2 
pulse  results.  Note  that  the  spectrograms  produced  using  a 
smaller  signal  slice  in  the  Fourier  transform  became  too 
noisy  for  interpretation,  and  larger  signal  slices  resulted  in 
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FIG.  7.  (Color  online)  Spectrogram  of  the  FDTD  predicted  held 
values  in  the  DNG  medium  measured  at  the  point  z  =  200Az'  for 
the  M  l  input  pulse,  where  Ac' =  50  jura, 
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FIG.  9.  (Color  online)  Spectrogram  of  the  FDTD  predicted  field 
values  in  the  DNG  medium  measured  at  die  point :  =  lOOAz  for  tine 
5-10-5  input  pulse,  where  Az  =  100  /tim. 


broad  areas  in  the  spectrogram  where  the  entire  signal  spec¬ 
trum  was  displayed.  The  use  of  larger  signal  slices  also 
shifted  the  estimates  for  the  arrival  times  of  various  frequen¬ 
cies.  The  Fourier  transforms  were  zero-padded  to  a  size  of 
4096  to  obtain  sufficient  resolution  (<1  GHz  sample  spac¬ 
ing)  near  the  center  frequency  (30  GHz)  of  the  signals. 

The  resulting  spectrograms  shown  in  Figs.  7-9  then  indi¬ 
cate.  respectively,  a  good  approximation  of  the  frequency 
components  present  at  the  location  z  =  0.0l  m  at  various 
points  in  time  for  the  1-1-1,  2-16-2,  and  5-10-5  cases.  The 
three-cycle,  1-1-1  pulse  spectrogram  in  Fig.  7  shows  a  nega¬ 
tive  slope.  The  higher  frequencies  propagate  faster,  hence 
they  arrive  sooner  than  the  lower  frequency  components. 
While  this  behavior  is  present  in  the  turn-on  and  turn-off 
portions  of  the  results  shown  in  Figs.  8  and  9,  it  is  clearer  in 


Fig.  8  since  the  2-16-2  pulse  has  faster  rates  of  change  that 
are  more  closely  related  to  the  1-1-1  case. 

A  plot  of  the  Fourier  transform  output  for  data  sections 
near  the  leading  edge,  central  portion,  and  trailing  edge  of 
the  FDTD  DNG  results  for  the  2-16-2  input  pulse  is  shown 
in  Fig.  10.  All  of  the  spectra  are  normalized  to  unity  for 
comparison  purposes.  This  figure  clearly  demonstrates  that 
the  leading  edge  of  the  pulse  has  significantly  higher  fre¬ 
quency  content,  the  central  portion  of  the  pulse  is  dominated 
by  the  CW  frequency,  and  the  trailing  edge  of  the  pulse  has 
significantly  lower  frequency  content.  Correlating  this  with 
Fig.  8,  one  then  sees  that  those  frequencies,  hence  the  lead¬ 
ing  front  of  the  pulse,  arrive  first  to  preserve  causality. 

To  further  investigate  this  dynamic  pulse  reshaping  be¬ 
tween  the  positive  index  of  refraction  (P1R)  high  frequency 
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FIG.  8.  (Color  online)  Spectrogram  of  the  FDTD  predicted  field 
values  in  the  DNG  medium  measured  at  the  point  z=200Az'  for 
the  2-16-2  input  pulse,  where  Ac'  =  50  /u.m. 
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FIG.  10.  (Color  online)  Normalized  frequency  composition  of 
the  (1)  leading  edge,  (2)  central  portion,  and  (3)  trailing  edge  of  the 
DNG  FDTD  result  measured  at  die  point  z  =  2()0Ac'  for  the  2-16-2 
input  pulse,  where  Az'  =  50  /x m. 
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FIG.  1 1.  (Color  online)  Spectrogram  of  the  difference  between 
the  dispersive  DNG  2-16-2  FDTD  result  shown  in  Fig.  8  and  the 
corresponding  approximate  analytical  solution.  Red/dark  region  on 
the  right  (bluc/dark  region  on  the  left)  represents  larger  FDTD  (ap¬ 
proximate  solution)  values. 

components  and  the  N1R  center  frequency  and  low  fre¬ 
quency  components,  the  differences  between  the  spectro¬ 
grams  of  the  FDTD  dispersive  DNG  results  and  the  approxi¬ 
mate  solutions  (23)  were  obtained.  These  spectrograms  were 
also  zero-padded  at  negative  times  to  clarify  the  frequency 
components  in  the  earliest  portions  of  the  signals.  The  2-16-2 
results  measured  at  z  =  0.01  m  are  shown  in  Fig.  II.  Al¬ 
though  the  FDTD  DNG  and  the  approximate  solution  signal 
fronts  arrive  at  the  same  time  and  their  CW  portions  are  in 
very  good  agreement,  the  difference  spectrograms  li  1  ter  out 
most  of  the  CW  behavior  and  reveal  the  importance  of  the 
apparent  time  delay  and  dispersive  properties  of  the  leading 
and  trailing  edges  of  the  FDTD  results.  One  finds  that  the 
CW  portions  of  the  FDTD  results  are  not  fully  germinated 
until  a  time  later  than  with  the  approximate  solution.  Simi¬ 
larly,  the  trailing  edge  of  the  FDTD  result  contains  informa¬ 
tion  that  is  delayed  relative  to  the  approximate  solution. 
These  delays  reinforce  similar  observations  made  with  re¬ 
gard  to  Figs.  5  and  6. 

By  aligning  the  CW  portions  of  the  FDTD  and  the  ap¬ 
proximate  solution  results  and  then  generating  the  difference 
spectrogram,  one  can  filter  out  all  of  the  CW  behavior  and 
highlight  the  remaining  low  and  high  frequency  behavior. 
The  spectrogram  for  the  2-16-2  input  pulse  results  measured 
at  £  =  0.01  m  is  shown  in  Fig.  12.  The  transition  from  the  off 
state  to  the  CW  state  in  the  FDTD  DNG  data  shows  the 
necessary  presence  of  higher  frequency  components  very 
early  in  time  to  reconstruct  the  pulse  front.  The  lower  fre¬ 
quencies  of  the  leading  edge  needed  to  complete  the  re¬ 
sponse  in  the  DNG  medium  arrive  later.  This  explains  the 
slower  turn -on  of  the  FDTD  results  in  comparison  to  the 
approximate  solution.  The  trailing  edge  of  the  DNG  FDTD 
result  does  contain  significantly  more  of  the  lower  frequency 
components  while  the  presence  of  higher  frequency  compo¬ 
nents  earlier  in  time  aids  the  transition  from  the  CW  portion 
to  the  off  state.  This  behavior  adds  delay  to  that  already 


FIG.  12.  (Color  online)  Spectrogram  of  the  difference  between 
the  dispersive  DNG  2-16-2  FDTD  result  shown  in  Fig.  8  and  the 
corresponding  approximate  analytical  solution  when  their  CW  por 
tions  are  aligned.  Red/dark  region  in  the  left-upper  and  right  lower 
portions  of  the  left  side  and  in  the  left-upper  and  right-lower  por¬ 
tions  of  the  right  side  of  the  figure  represents  larger  FDTD  values. 
Blue/dark  region  in  the  left-lower  and  right-upper  portions  of  the 
left  side  and  in  the  left-lower  and  right  upper  portions  of  the  right 
side  of  the  figure  represents  larger  approximate  solution  values. 

incurred  with  the  turn-on  portion  of  the  pulse  and  explains 
the  apparent  lag  in  the  turn-off  of  the  FDTD  results. 


V.  CONCLUSIONS 

The  one-dimensional  electromagnetic  problem  of  a  cur¬ 
rent  sheet  source  in  a  dispersive  DNG  medium  was  consid¬ 
ered.  A  lossy  Drude  model  of  the  DNG  medium  was  used. 
The  solution  was  generated  numerically  with  the  FDTD 
method.  The  analogous  problem  in  a  nondispersive  DNG 
medium  was  also  considered.  It  was  shown  that  the  solution 
to  this  problem  is  not  causaJ  in  agreement  with  similar  ob¬ 
servations  given  in  [5].  Approximate  solutions  that  combined 
causal  envelopes  and  the  NIR  properties  of  the  nondispersive 
DNG  medium  were  constructed:  they  were  shown  to  com¬ 
pare  well  to  the  FDTD  results  for  the  dispersive  DNG  me¬ 
dium.  It  was  thus  demonstrated  that  causal  results  do  indeed 
require  the  presence  of  dispersion  in  DNG  media.  Spectro¬ 
grams  of  the  FDTD  results  for  the  dispersive  DNG  medium 
and  of  the  differences  between  those  FDTD  results  and  die 
approximate  analytical  solution  revealed  that  the  highest  fre¬ 
quency  components,  which  experience  a  DPS  (PIR)  medium, 
arrive  first  and  are  responsible  for  generating  the  front  edge 
of  the  dispersive  DNG  pulse  results.  They  also  revealed  that 
the  trailing  edge  contained  significant  lower  frequency  com¬ 
ponents  that  see  a  DNG  medium  and  have  much  slower 
wave  speeds.  The  apparent  delays  in  the  germination  of  the 
CW  portion  of  the  pulse  and  in  its  termination  were  corre¬ 
lated  with  these  dispersive  effects.  Thus,  the  dynamic  reshap¬ 
ing  of  the  DNG  pulses  occurs  because  of  the  different  arrival 
times  for  die  higher  and  lower  frequency  components. 

It  was  also  demonstrated  that  the  CW  portions  of  the 
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pulses  do  obey  all  of  the  NIR  effects  expected  from  a  time- 
harmonic  analysis  in  a  nondispersive  DNG  medium.  Thus, 
CW  analyses  of  DNG  media  are  credible  as  long  as  very 
narrow  bandwidth  pulse  trains  are  considered  for  any  practi¬ 
cal  realizations.  This  has  been  the  case  in  all  of  the  experi¬ 
mental  results  reported  to  dale.  Moreover,  time  delays  for  the 
realization  of  the  NIR  effects  are  inherent  in  the  processes 
dictated  by  the  dispersive  nature  of  the  physics  governing 
them. 

Even  though  the  results  presented  here  were  derived  for 
the  broad  bandwidth  DNG  Drude  model,  it  is  anticipated  that 
the  dynamic  reshaping  due  to  dispersive  effects  will  also 
play  a  similar  role  in  current  metamaterial  realizations  of 
DNG  media.  This  is  particularly  true  of  the  planar  implemen¬ 
tations  [15-18]  that  realize  the  Drude  behavior.  It  should 
also  be  true  for  the  nonplanar  metamaterial  constructs  that 
involve  highly  resonant  elements  that  achieve  DNG  proper¬ 
ties  in  narrow  frequency  bands.  In  fact,  this  dynamic  reshap¬ 
ing  will  occur  at  a  faster  rate  since  changes  in  the  medium 
parameters  from  the  DPS  to  the  DNG  states  will  be  much 
mom  rapid.  Furthermore,  the  dynamic  reshaping  must  also 
occur  at  the  interfaces  between  DPS  and  DNG  media,  The 
delays  in  the  reflected  pulse  formation  observed  in  [12]  and 

[13]  occur  as  the  dispersive  effects  sort  out  the  frequency 
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components  for  the  wavefronts  and  the  appropriate  propaga¬ 
tion  directions  for  the  reflected  and  transmitted  waves.  Since 
most  of  the  current  nonplanar  DNG  metamaterial  experi¬ 
ments  are  of  the  slab  variety,  these  dynamic  reshaping  effects 
would  be  measurable  if  the  experimental  results  were  re¬ 
solved  in  time. 

It  is  interesting  to  note  that  the  DNG  results  presented 
here  can  be  related  to  Feynman’s  description  of  the  interac¬ 
tions  of  electrons  and  positrons  [28].  In  particular,  if  one 
views  the  causal  envelope  propagating  forward  in  time  as  an 
electron  and  the  NIR  CW  wave  propagating  backward  in 
time  as  a  positron,  then  the  positronlike  CW  wave  is  annihi¬ 
lated  by  the  electronlike  envelope  at  the  pulse  front.  The 
scattered  particle  is  the  resulting  wave  that  contains  the  NIR 
properties  but  is  causal  and  propagating  forward  in  time.  The 
vertex  interaction  is  smeared  in  time  by  the  dispersive  ef¬ 
fects;  it  funnels  the  highest  frequencies  into  the  pulse  front 
and  properly  distributes  the  higher  and  lower  frequencies  to 
provide  the  transition  from  the  off-state  to  the  CW  on -state. 
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The  nanotechnology  area  holds  much  promise  for  the  realization  of  ultra-small  and  ultra-fast 
devices  with  a  variety  of  interesting  applications.  The  optical  sector  is  an  immediate  beneficiary 
of  nanotechnologies.  Consider,  for  instance,  a  coated  nano-cylinder  or  nano-sphere  whose  outer 
radius  is  a  —  50  nm  .  At  an  excitation  wavelength  of  ^  =  500  nm ,  its  radius  is  only  a  =  !\  0 . 

One  could  use  arrays  of  such  nano-cylinders  or  nano-spheres  as  inclusions  in  a  substrate  to  realize 
optical  metamaterials.  On  the  other  hand,  because  of  their  subwavelength  sizes,  such  arrays 
could  be  used  for  several  nano-photonic  circuit  concepts. 

We  have  been  studying  analytically  and  with  the  finite  difference  time  domain  (FDTD)  approach, 
the  scattering  of  Gaussian  beams  from  sets  of  nano-cylinders  and  nano-spheres.  We  selected  the 
FDTD  approach  because  of  its  versatility  in  the  choices  and  configurations  of  materials  and 
structures  that  can  modeled.  This  versatility  makes  it  an  excellent  candidate  for  studying  the 
behavior  of  these  ultra-small  systems.  One  strong  design  goal  has  been  to  determine  whether  or 
not  an  array  of  plasmonic  nano-cylinders  could  be  used  as  an  effective  waveguide  to  localize 
efficiently  an  electromagnetic  field  to  subwavelength  dimensions. 

It  has  been  found  theoretically  and  experimentally  that  by  properly  taking  into  account  the  actual 
material  characteristics  at  optical  frequencies  and  by  properly  designing  the  nano-structures,  one 
can  achieve  enhanced  scattering.  In  particular,  the  permittivity  of  metals  such  as  gold  and  silver 
at  visible  (optical)  frequencies  can  be  described  effectively  by  a  lossy  Drude  (cold  plasma) 
model.  One  finds  that  the  real  part  of  the  permittivity  of  metals  is  in  fact  generally  negative  at 
those  frequencies.  On  the  other  hand,  it  is  readily  shown  from  canonical  plane  wave  scattering 
problems  that  the  total  scattering  cross-section  of  an  electrically  small  dielectric  cylinder  (sphere) 
(i.e.,  Rayleigh  scatterers  that  are  much  smaller  than  the  wavelength)  is  resonant  for  £  =  -l  .0£0 

(£  =  -2.0£o)  .  One  typically  realizes  these plasmon  resonance  effects  at  optical  frequencies  and 
at  nano-meter  length  scales  with  cylindrical  (spherical)  dielectric  scatterers  that  are  coated  with 
metals  whose  permittivity  is  €r  =  —  1 .0  (£r  =  —2.0)  at  the  wavelength  of  interest, 


We  have  used  a  finite  difference  time  domain  (TDTD)  simulator  to  study  the  enhanced  scattering 
from  plasmonic  cylinders  [1].  Here,  each  cylinder  is  defined  by  an  interior  free  space  (and  other) 
dielectric  cylinder  whose  radius  q  =  40 nm  that  is  coated  with  a  Drude  metal  shell  whose  outer 

radius  r2  =  50 nm  .  The  Drude  metal  coating  has  £r  =  -1 .0  at  =  500 nm ,  the  wavelength  of 

the  optical  field  under  consideration.  We  have  combined  together  an  array  of  these  cylinders  to 
form  a  nanometer  waveguide  environment  that  can  capture  a  significant  amount  of  power  from  an 
incident  optical  beam  and  can  then  transform  it  into  a  localized  field  with  a  highly  subwavelength 
spot  at  the  output  of  the  waveguide.  This  geometry  is  shown  in  Fig.  la.  A  parallel  polarized 
Gaussian  beam  (electric  field  in  the  plane  of  the  figure,  magnetic  field  orthogonal  to  it),  whose 


waist  is  2 Aq  =  1000  nm,  is  incident  normally  on  the  end  of  an  18  element  two  column  array  of 
these  plasmonic  nano-cylinders.  The  longitudinal  distance  between  the  centers  of  the  cylinders  is 
1 10  nm  giving  a  10 nm  separation  between  them  along  the  beam  axis.  The  entire  waveguide  is 

thus  slightly  longer  than  2 Aq.  The  transverse  distance  between  the  centers  of  the  cylinders  is 
120  nm  so  that  the  channel  formed  between  the  edges  of  the  cylinders  is  20  nm  .  The  distance 
between  the  outer  edges  of  the  cylinders  is  thus  220 nm  .  To  isolate  the  input  end  of  the 
waveguide  and  the  incident  beam  from  its  output  end,  the  cylinders  are  fed  through  a  slot  in  a 
1 0  nm  -thick  perfect  electric  conductor  (PEC)  ground  plane.  The  slot  width  is  1 70 nm  .  Detector 

arrays  (7/  elements,  each  element  being  6  nm  wide)  are  located  2 r2  =  100  nm  away  from  the 
ends  of  the  waveguide.  To  further  localize  the  field  at  the  output  end,  the  last  two  nano-cylinders 
have  10/7/n  slots  in  their  metallic  shells  [2,  3]. 

The  incident  beam  is  launched  from  a  total-field/scattered  field  boundary  into  the  FDTD 
simulation  region.  The  beam  scatters  from  and  excites  the  input  end  of  the  waveguide.  We  are 
investigating  the  waveguide  output  fields  as  functions  of  all  of  the  geometry  and  material  factors. 
Snapshots  in  time  of  the  magnetic  field  intensity  distribution  for  the  indicated  plasmonic  nano¬ 
cylinder  waveguide  configuration  are  shown  in  Fig.  1.  Hie  incident  field  is  clearly  coupled  to  the 
waveguide  and  the  transmitted  field  is  basically  confined  to  the  regions  between  the  nano¬ 
cylinders.  Each  nano-cylinder  is  also  seen  to  be  strongly  coupled  to  its  nearest  neighbors. 
Transverse  cross-sections  of  the  magnetic  field  intensity  at  various  locations  along  the  beam  axis 
are  shown  in  Fig.  2.  The  cross-sections  before  and  just  into  the  waveguide  are  given  in  Fig.  2a; 
those  just  after  the  output  of  the  waveguide  are  given  in  Fig.  2b.  The  full  width  at  half  maximum 
of  the  intensity  of  the  input  beam  is  235  nm  .  [t  is  found  that  the  intensity  distribution  across  the 
front  detector  array  has  a  145  nm  full  width  at  half  maximum  while  it  has  a  1 20  nm  full  width  at 
half  maximum  across  the  back  detector  array.  The  output  beam  is  seen  to  be  highly  localized 
near  the  waveguide.  It  has  a  80  nm  full  width  at  half  maximum  across  the  ouput  face  of  the 
waveguide.  This  spot  size  is  6.25  time  smaller  than  the  wavelength.  From  the  time  histories 
obtained  at  the  center  of  the  front  and  back  detectors  it  is  found  that  the  amplitude  of  the  magnetic 
field  at  the  back  detector  array  is  about  20%  of  its  value  at  the  front  detector  array  and  that  it 
takes  only  about  15  cycles  for  this  level  to  be  reached.  The  plasmonic  nano-cylinder  waveguide 
array  has  thus  localized  the  incident  field  with  a  reasonable  efficiency  to  subwavelength  spot 
sizes.  These  localized  fields  could  be  used  for  many  nanoscale  optical  probe  and  photonic  circuit 
applications. 
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Fig.  1 .  Snapshots  in  time  of  the  magnetic  field  intensity:  a)  t  =  0  At,  b)  t  =  150  At,  c)  t  = 
1200  At,  and  d)  t  =  3200  At.  The  field  translated  along  the  plasmonic  nano-cylinder 
waveguide  produces  a  localized  distribution  at  the  output  face. 


FDTD  cell  number 


Fig.  2.  Transverse  cross-sections  of  the  magnetic  field  intensity  distribution  at  the  (a) 
input  and  (b)  output  ends  of  the  plasmonic  nano-cylinder  waveguide. 


